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Executive Summary 

Simulation of turbulent mixing and chemical processes in the near-field plume 
and plume- vortex regimes has been successfully carried out recently using a reduced gas 
phase kinetics mechanism which substantially decreased the computational cost. A 
detailed mechanism including gas phase HOx, NOx, and SOx chemistry between the 
aircraft exhaust and the ambient air in near-field aircraft plumes is compiled. A reduced 
mechanism capturing the major chemical pathways is developed. Predictions by the 
reduced mechanism are found to be in good agreement with those by the detailed 
mechanism. With the reduced chemistry, the computer CPU time is saved by a factor of 
more than 3.5 for the near-field plume modeling. Distributions of major chemical species 
are obtained and analyzed. The computed sensitivities of major species with respect to 
reaction step are deduced for identification of the dominant gas phase kinetic reaction 
pathways in the jet plume. 

Both the near field plume and the plume-vortex regimes were investigated using 
advanced mixing models. In the near field, a stand -alone mixing model was used to 
investigate the impact of turbulent mixing on the micro- and macro-scale mixing 
processes using a reduced reaction kinetics model. The plume-vortex regime was 
simulated using a large-eddy simulation model. Vortex plume behind Boeing 737 and 
747 aircraft was simulated along with relevant kinetics. Many features of the computed 
flow field show reasonable agreement with data. The entrainment of the engine plumes 
into the wing tip vortices and also the partial detrainment of the plume were numerically 
captured. The impact of fluid mechanics on the chemical processes was also studied. 
Results show that there are significant differences between spatial and temporal 
simulations especially in the predicted S0 3 concentrations. This has important 
implications for the prediction of sulfuric acid aerosols in the wake and many partly 
explain the discrepancy between past numerical studies (that employed parabolic or 
temporal approximations) and the measured data. 

Finally to address the major uncertainty in the near-field plume modeling related 
to the plume processing of sulfur compounds and advanced model was developed to 
evaluate its impact on the chemical processes in the near wake. A comprehensive aerosol 
model is developed and it is coupled with chemical kinetics and the axisymmetric 
turbulent jet flow models. The integrated model is used to simulate microphysical 
processes in the near-field jet plume, including sulfuric acid and water binary 
homogeneous nucleation, coagulation, non-equilibrium heteromolecular condensation, 
and sulfur-induced soot activation. The formation and evolution of aerosols are computed 
and analyzed. The computed results show that a large number of ultra-Fine (0.3-0.6 nm 
in radius) volatile H 2 S0 4 - H 2 0 embryos are generated in the near-field plume. These 
embryos further grow in size by self coagulation and condensation. Soot particles can be 
activated by both heterogeneous nucleation and scavenging of H 2 SO A - H 2 0 aerosols. 
These activated soot particles can serve as water condensation nuclei for contrail 
formation. Conditions under which ice contrails can form behind aircrafts are studied. 
The sensitivities of the threshold temperature for contrail formation with respect to 
aircraft propulsion efficiency, relative humidity, and ambient pressure are evaluated. 
The computed aerosol properties for different extent of fuel sulfur conversion to S(VI) 
(S0 3 and H 2 S0 4 ) in engine are examined and the results are found to be sensitive to this 
conversion fraction. 



1. PROBLEM DESCRIPTION 

The projected increase in future aviation especially in heavy traffic flight 
corridors has the potential for causing a major upper atmospheric perturbation. Recent 
research has focused on the potential impact on the upper troposphere and on the lower 
stratosphere due to exhaust emission from the proposed fleet of supersonic high-speed 
civil transport (HSCT) and/or advanced subsonic transport (AST) operating in the 
transatlantic and Pacific rim flight corridors (e.g., Brasseur et al., 1990, Beck et al., 1992; 
Stolarski and Wesoky, 1995; Friedl, 1997). Emissions from individual aircraft are 
initially highly localized in the engine plume wake and rapidly disperses in the far wake. 
However, the chemical processes in the hot species-rich plume may lead to changeson 
the global scale if it causes marked changes in the far field exhaust composition, inis 
concern, especially as related to the potential deployment of HSCT and/or AST fleets, has 
led to many investigations focusing on the near field (characterized as the exhaust plume 
and plume-wing-tip vortex interaction regimes) chemical processes. Gas phase chemical 
processes due to NO x , HO x and SOx species and heterogeneous kinetics due to exhaust 
aerosols and soot particles in the plume are being studied (Fahey et al., 1995a, 1995b; 
Brown et al., 1996a, b; Weisenstein et al., 1993). Both measurements and numerical 
studies in the near field indicate potentially significant impact on the atmosphere 
depending upon the significance of both homogeneous and heterogeneous kinetics. In 
particular, conversion of emitted NO x to inactive forms of nitrogen, plume processing of 
sulfur and formation of ice in the plume have been the focus of many recent 
investigations (Friedl, 1997). 

The net effect of NO x emission on ozone is known to be reduced if some of it is 
converted to its less reactive form (via the hydrolysis of N 2 O s to HN0 3 ). Near-field 
modeling studies (Miake-Lye et al., 1993, 1994; Quackenbush et al., 1993) and 
measurements (Fahey et al., 1995a) indicate that relatively small amounts of NO y is 
deposited in the form of nitrous and nitric acid with the bulk remaining in the form of 
NO x . Thus, both data and measurements suggest that there is only a minimal conversion 
of emitted NO x into other species in the near field plume. On the other hand, the sulfur 
compounds in the near field plume are known to be significantly affected by the chemical 
processes. Recent calculations (Miake-Lye et al., 1994, Brown et al., 1996a, d) and 
experimental observations (Fahey et al., 1995b; Karcher and Fahey, 1997) suggest that 
many small aerosol particles are formed when H 2 S0 4 (v) undergoes binary homogeneous 
nucleation with water. The number of these particles is estimated to be orders of 
magnitude larger than the emitted soot particles. Calculations (Miake-Lye et al., 1994) 
suggests that the nucleation process is much more likely than heterogeneous condensation 
of H 2 S0 4 on soot particles. The resulting increase in overall aerosol surface area in the 
plume has the potential for changing the exhaust composition in the far field wake and 
thereby, impacting the global dynamics. 

The major uncertainty is in the modeling of soot and aerosol particles in the 
plume. For example, measurements (Fahey et al., 1995b) suggests that around 10% of the 
fuel sulfur is converted to sulfates whereas models (Brown et al., 1996a, c) predict only 
around 1-2 %. This discrepancy has been the focus of recent studies (e.g., Brown et al., 
1996a-d; Schumann et al., 1996; Karcher and Fahey, 1997) and is of fundamental interest 
since the formation of large number of aerosol particles in the plume may lead to 
significant global effects in the long term. For reliable estimates, physical processes that 
involve gas-particle condensation, homogeneous binary nucleation of H 2 SO A — H 2 0 
aerosols, condensation of H 2 SO A - H 2 0 on soot particles, aerosol coagulation and soot 



agglomeration, etc., all must be modeled. Although most of these processes have been 
modeled in the past, either separately or together, the coupling of these processes along 
with the gas phase kinetics in a fluid dynamic solver has been difficult primarily due to 
the computational cost. Both condensation aerosols and coagulation of aerosols require 
that the number density and size distribution of the particles be determined as a function 
of space and time (due to their dependence on the local composition and thermodynamic 
state). The computational cost to track all the aerosol sizes is excessive even for 2D 
steady state calculations (although, recently Brown et al., 1996a, d have begun doing 
some calculations). Evaluating aerosol dynamics in an unsteady reacting environment by 
using LES is impossible at present. Unsteady mixing effects on aerosol dynamics may be 
important since condensation and coagulation processes could be modified if aerosol or 
soot particles are entrained in local vortical structures containing large concentrations of 
key species (e.g., OH,SO^,H 1 0). Reliable predictions of plume processing of sulfur may 
require taking all these issues into account. 

The exorbitant cost of computing the details of aerosol dynamics appears to be the 
main stumbling block in the study of heterogeneous kinetics in the near-field plume. 
However, if the models used for condensation and coagulation can be reformulated so 
that a significant reduction in computational cost is achieved without causing a major loss 
of accuracy, then it would be feasible to address these processes within the context of 
even an unsteady simulation of the near-field. The advantages of a computationally 
efficient model for aerosol dynamics are obvious. In the present study, both 
homogeneous and heterogeneous kinetics modeling issues have been addressed with the 
view for application in large-scale 3D simulation methodology such as large-eddy 
simulations. 

This report summarizes the results obtained in this study. Most of the details are given in 
the papers attached as Appendices. Here we briefly summarize the pertinent results. 


2. SUMMARY OF RESULTS 

Many past studies (Miake-Lye et al., 1993; 1994; Quackenbush et al., 1996; 
Brown et al., 1996a) of the near field plume have employed a variant of the standard 
plume code (SPF-II) called UNIWAKE. This code is a space-marching, parabolized 
solver and employs the standard k-e turbulence model to obtain the steady state field. 
The computational efficiency of this code is achieved by using a parabolic (2D) approach 
and by modeling all chemical and fluid dynamic processes in the steady state limit. As a 
result, physical time dependent processes that govern mixing and combustion processes 
such as molecular diffusion, large-scale entrainment and small-scale turbulent mixing are 
ignored (or modeled using eddy viscosity/diffusivity models). Thus, the importance (if 
any) of the unsteady processes on the chemical processing in the plume could not be 
addressed using this type of code. 

To determine the impact of the time-dependent processes on plume dynamics, a 
study was undertaken under this NASA AEAP program. Since computations that resolve 
the unsteady processes are extremely expensive, new cost-effective methodologies had to 
be developed and incorporated for this study. 

The project had four primary objectives: 

(1) to develop and validate reduced chemical kinetics models for near- field studies 

(2) to evaluate micro- and macro-scale mixing effects in the plume regime 

(3) to investigate the effects of mixing on chemical processes in the plume-vortex regime 

(4) to investigate reduced (computationally efficient) models for heterogeneous kinetics 



We have achieved all of the above noted objectives and have obtained new insight 
into the fluid dynamic and chemical processes occurring in the near field. The 
investigation to develop reduced heterogeneous models is currently underway. However, 
the complexity of the heterogeneous processes that must be included, the lack of detailed 
data for comparison and the computational cost of doing baseline (i.e., conventional 
discrete-nodal/sectional method) simulations has delayed the development effort. 
Currently, we are including aerosol models into our codes and we expect to get some 
results before the end of the current project. However, discrete-nodal/sectional method is 
computationally too expensive for general purpose studies. To develop computationally 
efficient aerosol dynamics models, we have identified two methodologies that show 
potential for application within fluid dynamic solvers. However, the various issues that 
must be addressed cannot be completed under the current program. In the following, we 
briefly summarize the key results of the current research. More details are given in 
Menon and Chen (1995), Menon and Wu (1998) Wu and Menon (1998a, 1998b), Wang 
and Chen (1997). 

2. 1 Reduced Kinetics Development 

Reduced reaction kinetics are modified chemical reaction models that mimic 
accurately all the critical features of a full (or detailed) reaction set. An obvious 
advantage of reduced kinetics is that it provides a computationally affordable set of key 
species that can be used in unsteady models and in large-eddy simulations (LES). Gas- 
phase mechanisms appropriate for the engine exhaust plume in the upper troposphere and 
lower stratosphere were compiled from literature. A 25-species, 47 reaction set used 
earlier for the HSCT (Miake-Lye et al., 1993; Brown et al„ 1996) was reduced to a 8- 
species, 16 reaction set. Note that this approach does not just neglect the other reactions 
and species. Rather, it modifies the reaction rates in the reduced set using sensitivity 
analysis so that the resolved species production/destruction time scales reflect the impact 
of the neglected species (which are assumed to be in quasi-equilibrium). 

The dominant oxidation kinetic pathways are identified through sensitivity analysis of 
major species with respect to gas phase reaction steps. The spatial distributions of major 
chemical species in the near-field jet plume are computed and analyzed. The reduced 
mechanism is developed which captures the major chemical pathways in the near-field jet 
plume. Results obtained with the reduced chemistry show good agreement with those 
predicted by the detailed chemistry. However, using the reduced chemistry has cut down 
computer CPU time by a factor of more than 3.5 for the near-field simulation. Results 
also show that the reduced mechanism captures the major pathways responsible for ozone 
destruction very accurately and also predicts accurately all the species that were included 
in the reduced set (Menon and Chen, 1995; Wang and Chen, 1997). 

2. 2 Turbulent Mixing and Chemical Processes in the Plume Regime 

The potential impact of engine exhaust gas on ambient ozone has been studied 
using several numerical models with differ in their treatments of the mixing process 
between the exhaust gas and the ambient air. These models include Perfectly-Stirred 
Reactor (PSR), transient Well-Mixed Reactor (WMR), Partially-Stirred-Reactor (PaSR), 
joint scalar Probability Density Function (PDF) and the Linear-Eddy Mixing (LEM) 
models. The last two methods, i.e., PDF and LEM (developed during this study for near- 
field plume studies) are time dependent simulation methods that incorporate effects of 
both micro- and macro-scale mixing. The PDF model is similar to the standard plume 
model in terms of the fluid dynamics but it incorporates a Monte-Carlo stochastic 



simulation for the joint scalar pdf. This allows proper incorporation of the effects of 
turbulence on the chemical process. However, only the LEM approach includes micro- 
scale mixing using physically consistent turbulent scaling rules and incorporates, in 
addition the effects of molecular diffusion. The results for an exhaust plume of a typical 
HSCT aircraft (Menon and Wu, 1998) showed that compared to predictions based on 
models that do not include small-scale mixing and/or molecular diffusion effects, 30% 
less Oj and 15% less NO x are depleted in the near field of the plume. This suggests that 
local unmixedness can have an inhibiting effect on ozone depletion in the near-field 
plume. Inclusion of the formation of condensed nitric acid on soot particles showed that 
15% of the available N0 X is converted into its inactive form but has negligible effect on 
O3 concentration. Competition between the gas phase H 2 S0 4 formation and the H 2 S0 4 
condensation on soot particles for the available S0 3 in the plume was found to depend 
partly on the mixing process. Results (see Menon and Wu, 1998) show that for the single 
plume case, the LEM results are quite similar to the earlier models that ignored or 
modeled in an ad hoc manner, the mixing processes. This suggests that the macro (large) 
scale entrainment process (included in various models) dominate the mixing process. 
However, as shown in Wang and Chen (1997) and Menon and Wu (1997) there are 
significant differences in the local details of the radical (0H,S0 2 ,S0 3 ) distribution that is 
likely to impact the heterogeneous kinetics related to sulfuric acid. Effects of unsteady 
mixing on the aerosol and soot related chemical processes remains unresolved at this 
time. Future study will address this particular issue. 


2.3 LES of Reacting Plume-Vortex Interactions 

The development of the reduced reaction mechanism provided a computationally 
affordable model for LES. LES is probably the only method available to study unsteady 
mixing and chemical processes during plume- vortex interaction. In LES, all scales larger 
than the grid are resolved using a time- and space-accurate scheme and only the small 
scales are modeled. We have developed a new localized dynamic subgnd model for the 
subgrid kinetic energy to close the LES equations. Using this approach with a 15-species 
reduced gas phase reaction mechanism (and the heterogeneous kinetics involving nitric 
and sulfuric acid condensation on soot particles) a series of calculations were earned out 
to determine (a) the dynamics of fluid mixing and entrainment of the engine exhaust into 
the wing tip vortex, (b) the impact of mixing on the kinetics in the plume-vortex stage, 
and, (c) the capability of temporal versus spatial LES to capture the entrainment and the 

chemical processes. 

The 3D LES was carried out using both temporal and spatial methods. In temporal 
simulation, periodic boundary conditions are used in the streamwise direction and the 
simulation captures a region of plume-vortex evolving in time which corresponds 
spatially (via a velocity transformation) to a location that moves (with increasing time) m 
the downstream direction. An advantage of this approach is that simulations can be used 
to study the long time (far downstream) evolution of the plume-vortex interaction. This 
approach is always used in the far field wake break up regimes (e.g., Lewellen and 
Lewellen, 1997) and in global modeling. However, temporal simulations cannot capture 
the inflow-outflow effects which may be important in the initial stages of the interaction 
between the engine exhaust plume and the wingtip vortices. Of particular interest here is 
the process of entrainment and detrainment (whereby portion of the plume is not 
entrained into the vortex, as observed in flight). To investigate these issues, two 
simulations using identical resolution and conditions were carried out using temporal and 
spatial methods. For the latter method, proper inflow turbulence and non-reflective 
outflow conditions were implemented. Reduced kinetics including the heterogeneous 



processes were studied for an initial period (approximately 6-8 seconds) of the plume- 
vortex interaction phase. Although these simulations are quite expensive, using an 
optimized massively parallel code, a spatial LES using 1.2 million grid points and with 
15 reacting species takes approximately 2000 single-processor hours on a SGI Origin 
2000 to simulate 5 flow-through times. This CPU requirement is not considered 
exorbitant considering the availability of fast parallel systems. 

Results show that the entrainment of the plume into the vortex in both the Boeing 
737 and 747 cases are qualitatively similar to the LIDAR measurements. Results also 
showed that for the 747 case, a portion of the engine plume is not entrained into the 
vortex and if fact, appears to split away (a form of detrainment similar to that observed 
earlier in flight). To evaluate the effect of modeling strategy on the mixing process 
temporal (as is typically used by most researchers) and true spatial evaluation were 
simulated under otherwise identical conditions. Although both temporal and spatial LES 
look similar there are significant differences in the entrainment process and in the actual 
distribution of typical species. The breakup of the plumes into two pieces (one that is 
entrained and another that is left behind) is seen to occur more readily in the spatial case. 
More importantly, estimated of the mass of NO and S0 3 entrained into the vortex core 
show that the temporal simulation overpredicts the entrainment of NO by nearly 25% by 
4 seconds while under predicting entrained S0 3 by nearly 25%. Since S0 3 is formed due 
to gas phase reactions and subsequently plays a major role in the formation of sulfuric 
acid aerosols as well as in the formation of condensed sulfuric acid on soot particles, this 
result has some significant implication when dealing with heterogeneous effects. 

In summary, unsteady mixing and gas phase kinetics have been incorporated into 
two different types of fluid dynamic solvers in the plume regime: the PDF and LEM 
methods. Inclusion of aerosol dynamics into these codes using the discrete-sectional 
method is currently underway. The plume-vortex interaction regime has been addressed 
using dynamic subgrid model based LES. Inclusion of gas phase kinetics and simple 
heterogeneous processes have been successfully carried out. Details of the mixing and 
chemical processes have been extracted from these simulations and results point to the 
importance of dealing with inflow/outflow processes when simulating the initial process 
of plume-vortex interaction. The computational cost has been made reasonable using an 
optimized parallel version of the 3D code. 


2.4 Aerosol and Soot Dynamics Modeling in the Plume Regime 

A detailed chemical mechanism between aircraft exhaust and ambient air in near- 
field jet plumes was earlier developed and compiled. To further extend the chemistry 
capability, a comprehensive aerosol model is developed to simulate micro physical 
processes including sulfuric acid and water binary homogeneous nucleation, aerosol 
coagulation, non-equilibrium heteromolecular condensation, and sulfur-induced soot 
activation as water condensation nuclei. These micro physical processes are coupled with 
gas-phase chemistry and used to study axisymmetric turbulent jet flow. 

The modeling results indicate that a large number of subnanometer-sized volatile 
H 2 SC >4 - H 2 0 embryos are generated in the near-field plume. These embryos grow in 
size by condensation in the early jet regime and by self coagulation later. Soot particles 
can be activated into water condensation nuclei in the jet plume for possible contrail 
formation through both heterogeneous nucleation and scavenging of volatile 

H 2 S0 4 - H 2 0 aerosols. 



The gas phase conversion of S0 2 to H 2 S0 4 is controlled by OH driven oxidation 
and it is limited to less than 1%. However, the amount of H 2 S0 4 produced will increase 
via fast reaction of S0 3 with H 2 0 directly if S0 3 emission due to fuel sulfur oxidation 
to S0 3 in the engine is considered. Higher levels of H 2 S0 4 in turn result in more sulfate 
aerosol embryos formed in the early jet stage, and also enhance condensation of 
H 2 S0 4 and H 2 0 onto the sulfate aerosols in this regime. 

Consequently, S0 3 emission from the engine nozzle exit increases the 
cumulative number density and surface area density of volatile sulfate aerosols. S0 3 
emission also promotes chemical activation of soot particles with enhanced acid coating. 
This suggests the need for measurements of the partitioning of sulfur species ( S0 2 and 

S0 3 ) at the exhaust exit. 

The threshold condition for ice contrail formation behind an aircraft depends on 
the ambient pressure, the ambient relative humidity, the overall propulsion efficiency of 
the aircraft and the fuel properties such as the fuel hydrogen mass fraction and the 
specific combustion heat. Higher values of ambient relative humidity or pressure 
facilitate contrail formation. The threshold temperature is found sensitive to the overall 
propulsion efficiency of the aircraft. The prediction of particle and species distribution in 
the near-field plume can serve as input for the follow-on larger-scale model. The present 
work also provides an estimate of the aerosol surface area density for subsequent 
evaluation of heterogeneous oxidation mechanisms in the wake. 
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ABSTRACT 

Turbulent mixing and chemical reactions in the near field of an engine exhaust jet plume have been investigated 
using a mixing model that explicitly incorporates both large- and small-scale turbulent mixing and the molecular 
diffusion effects. A reduced reaction mechanism that accurately reproduces results of a detailed mechanism for 
the lower stratosphere was used to simulate the exhaust plume dynamics of a typical high-speed civil transport 
aircraft. This study shows that, compared to predictions based on models that do not include small-scale mixing 
and/or molecular diffusion effects. 30% less O, and 15% less NO, are depleted in the near field of the plume 
This suggests that the lack of local mixing can have an inhibiting effect on ozone depletion in the near-field 
plume. Inclusion of heterogeneous kinetics involving the formation of niiric and sulfuric acid due to water 
condensation on soot particles showed that 15% of the available NO, is converted into its inactive form but has 
a negligible effect on O, concentration. This method also provides appropriate conditions for the plume-vortex 
interaction stage where more complex chemistry, including heterogeneous kinetics, is likely to take place. 


1. Introduction 

The projected increase in future aviation, especially 
in heavily trafficked flight corridors, has the potential 
for causing a major upper-atmospheric perturbation. Re- 
cent research has focused on the potential impact on the 
upper troposphere and on the lower stratosphere due to 
exhaust emissions from the proposed fleets of super- 
sonic high-speed civil transport (HSCT) and/or ad- 
vanced subsonic transport (AST) operating in the trans- 
atlantic and Pacific rim flight corridors (e.g., Brasseur 
et al. 1990; Beck et al. 1992; Stolarski and Wesoky 
1995; Friedl 1997). Emissions from individual aircraft 
are initially highly localized in the engine plume wake 
and rapidly dispersed in the far wake. However, the 
chemical processes in the hot species-rich plume may 
lead to changes on a global scale if this causes marked 
changes in the far-field exhaust composition. This con- 
cern, especially as related to the potential deployment 
of HSCT and/or AST fleets, has led to many investi- 
gations focusing on the near-field (characterized as the 
exhaust plume and plume-wing tip vortex interaction 
regimes) chemical processes. Gas phase chemical pro- 
cesses due to NO, (NO, N0 2 ), HO r (OH, H0 2 , H 2 0 2 ), 
and SO, (SO, S0 2 , S0 3 , HS0 3 , H 2 S0 3 ) species and 
heterogeneous kinetics due to exhaust aerosols and soot 
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particles in the plume are being studied (Fahey et al. 
1995a,b; Brown et al. 1996a,b; Weisenstein et al. 1996), 
Both measurements and numerical studies in the near 
field indicate a potentially significant impact on the at- 
mosphere, depending upon the significance of both ho- 
mogeneous and heterogeneous kinetics. In particular, the 
conversion of emitted NO, to inactive forms of nitrogen 
and the impact of plume processing of sulfur due to the 
formation of a large number of aerosol particles in the 
plume have been the focus of many recent investigations 
(Friedl 1997; Stolarski et al. 1995; WMO 1995). 

The net effect of NO, emissions on ozone is known 
to be reduced if some of it is converted to its less reactive 
form (via the hydrolysis of N 2 0 5 to HN0 3 ). Near-field 
modeling studies (Miake-Lye et al. 1993, 1994; Quack- 
enbush et al. 1993) and measurements (Fahey et al. 
1995a) indicate that relatively small amounts of NO, 
are deposited in the form of nitrous and nitric acid, with 
the bulk remaining in the form of NO,. Thus, both data 
and measurements suggest that there is only a minimal 
conversion of emitted NO , into other species in the near- 
field plume. On the other hand, the sulfur compounds 
in the near-field plume are known to be significantly 
affected by the chemical processes. Recent calculations 
(Brown et al. 1996a, b) and experimental observations 
(Fahey et al. 1995b; Karcher and Fahey 1997) suggest 
that many small particles are formed in the plume via 
binary homogeneous nucleation of sulfuric acid with 
water. This increase in the overall aerosol surface area 
in the plume has the potential for changing the exhaust 
composition in the far-field wake. 
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Many past studies (Miake-Lye et al. 1993, 1994; 
Brown et al. 1996a) of the near-field plume have em- 
ployed the standard plume code, SPF-II (Dash et al. 
1990). SPF-II is a space-marching, ax i symmetric solver 
and employs the standard k - e turbulence model to 
obtain the steady-state field. The computational effi- 
ciency of this code is achieved by using an axisymmetric 
approach and by modeling all chemical and fluid dy- 
namic processes in the steady-state limit. Farther down- 
stream, the evolution of the plume and its interaction 
with the wingtip vortices are modeled using codes such 
as UNIWAKE (Quackenbush et al. 1993, 1996) which 
uses the output from the SPF-II calculation of the plume 
as inflow conditions. UNIWAKE is a 3D unsteady, par- 
abolized Navier-Stokes solver that incorporates the ef- 
fect of turbulence using the second-order Reynolds 
stress closure. Detailed gas phase kinetics have been 
studied using this code. However, since all turbulent 
scales are modeled by the closure (in both SPF-II and 
UNIWAKE codes), many physical time-dependent pro- 
cesses that govern mixing and combustion processes 
(such as molecular diffusion and small-scale turbulent 
mixing) are ignored (or modeled) in these codes. To 
include the time-dependent physics, methods such as 
large-eddy simulations (LES) can be employed. Al- 
though LES has been successfully used in the far field 
with limited (i.e„ one-step chemistry) or no chemistry 
modeling (e.g., Sykes et al. 1992; Lewellen and Lew- 
ellen 1996), this approach is not practical for the near- 
field studies due to the computational cost (orders of 
magnitude larger than the parabolized methods) since a 
large number of species must be simulated along with 
the fluid dynamics. It would be computationally feasible 
if a nonreacting LES flow field was used as a repre- 
sentative fluid dynamic snapshot and the chemical re- 
actions were simulated assuming no temporal depen- 
dence between the fluid motion and the chemical pro- 
cesses. However, this approach would not resolve the 
importance (if any) of the unsteady mixing processes 
on the chemical processing in the plume. In addition, 
decoupled LES would not allow for proper transport of 
the various species and would result in potentially in- 
correct estimation of the chemical effect. Due to these 
factors, it appears that no study has been reported on 
the effects of unsteady turbulent mixing in the near-field 
plume dynamics to date. 

This paper describes an innovative simulation meth- 
odology that allows proper characterization of the effect 
of small-scale mixing processes on the chemical processes. 
Inclusion of microscale mixing (missing in the earlier 
plume studies) should result in a physically consistent sim- 
ulation model and should provide a more accurate estimate 
of the radical (e.g., OH, NO,, CO, SO r ) reactions that are 
highly dependent upon the mixing process. The model 
described in this paper allows incorporation of detailed 
kinetics (if warranted) without making the computational 
cost unduly exorbitant. The key innovative feature of this 
approach is that the microscale mixing processes (diffusion 


and turbulent stirring) are explicitly simulated at their char- 
acteristic timescales. Thus, turbulence-chemistry interac- 
tion evolves in a consistent manner and its effects on the 
scalar fields are obtained by statistically averaging over 
the period of simulation. This is in contrast to steady-state 
methods where all time -dependent processes are modeled. 

The domain simulated in this study extends from the 
engine nozzle exit to a point where the plume begins to 
interact with the wingtip vortex. This domain spans a 
distance of about 1 km and simulates the plume for about 
1-2 s. Although this spatial and temporal extent is very 
small relative to the overall scales for plume evolution, 
this near-field region has been identified as the region 
where significant chemical processing occurs (Brown et 
al. 1996a). The engine exhaust conditions are chosen for 
a representative HSCT aircraft flying at an altitude of 18 
km along with appropriate ambient conditions (Brown et 
al. 1996a; Miake-Lye et al. 1993). Both homogeneous 
and heterogeneous kinetics in the near-field plume are 
investigated. Gas phase kinetics involving HO ( , NO ( , and 
SO, reactions are simulated using a reduced reaction 
mechanism that significantly reduces the computational 
cost without losing any accuracy in the prediction of the 
important species. Condensation of water on emitted soot 
particles and heterogeneous kinetics resulting in con- 
densed HNO, and H : S0 4 formation on the surface of 
condensed solid water are investigated in this study. More 
detailed heterogeneous kinetics related to soot adsorption, 
homogeneous binary nucleation of sulfur particles, vol- 
atile H : 0-H : S0 4 condensation on aerosol particles, and 
aerosol coagulation (Brown et al. 1996a) are deferred for 
a later study. 

The unsteady mixing model demonstrated here is ex- 
pected to provide a more accrurate estimate of chemical 
processing in the species-rich plume when compared to 
SPF-II. Thus, the outflow (time-averaged data) from this 
model can provide a more accurate inflow condition for 
the regime where the plume interacts with the wingtip 
vortices. UNIWAKE calculations or LES of this regime 
can employ this data. Currently, LES in the plume-vortex 
regime is using this initialization approach (Menon and 
Wu 1997). 

2* Turbulent mixing models 

In this study, the predictions using the new mixing mod- 
el are compared to the predictions made using other well- 
known mixing models. These models are described below. 

a. Perfectly stirred reactorfaell-mixed reactor 
(PSR/WMR) 

This model is also sometimes called the “continuously 
stirred" reactor or “box" model (e.g., Kee et al. 1989; 
Karcher 1995; Danilin et al. 1994) and is the most com- 
monly used model (in one variation or other) in aircraft 
plume studies. In this model, there is a continuous mass 
flow in and out of a constant volume reactor. Within this 
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reactor, the scalar fields are assumed perfectly mixed so 
that the rate of chemical processes is controlled by the 
kinetics of the reaction mechanism and not by the mixing 
process. The key parameter for the calculation is the res- 
idence time in the reactor that determines the final reacted 
state of the initially specified scalar fields. These models 
cannot account for turbulent mixing by both the large and 
the small scales and. therefore, cannot be used to study 
turbulence effects. On the other hand, due to its compu- 
tational efficiency, PSR/WMR codes can be used to in- 
vestigate detailed mulfispecies reaction mechanisms very 
quickly and to evaluate reduced kinetics. The mathematical 
and computational aspects of this modeling approach are 
given in numerous past studies (e.g., Kee et al. 1989; 
Danilin et al. 1994) and, therefore, are omitted for brevity. 

b. Probability density function (PDF) simulation 

method 

In this approach, the joint scalar PDF for the scalar 
fields is solved by a transport equation using the Monte 
Carlo technique. The mean flow is modeled by the Reyn- 
olds stress closure with a downstream marching algo- 
rithm (Wang and Chen 1997) as in the standard plume 
codes (e.g., SPF-II). Although turbulence effects are 
explicitly included in the direct estimation (i.e., without 
any modeling) of the chemical source term, the turbulent 
mixing is included using an ad hoc model and molecular 
diffusion is ignored. This approach was used extensively 
earlier (Menon and Chen 1995; Wang and Chen 1997) 
to evaluate chemical kinetics. In this paper, we will 
compare predictions using the new mixing model 
against the PDF predictions. 

c. The linear-eddy model (LEM) for scalar mixing 

This model explicitly treats three physical processes: 
1 ) large-scale turbulent entrainment, 2 ) small-scale tur- 
bulent mixing, and 3) reaction-diffusion at the small 
scales as independent processes that occur at their char- 
acteristic times. In the following, a brief overview of 
these processes are given. Subsequently, in section 2 d, 
more details pertaining to the actual implementation for 
exhaust plume simulations are provided to identify the 
new and innovative features of this modeling approach. 

In LEM, the scalar reaction-diffusion equations are 
solved on a ID domain (and hence, called linear eddy) 
using an explicit finite-difference scheme on a high- 
resolution grid so that all the length scales (e.g., from 
the integral scale L to the Kolmogorov scale 17 ) are 
resolved. Thus, the reaction-diffusion evolution is con- 
sidered exact and requires no closure. Although this 
evolution is simulated in ID, this domain is not a phys- 
ically fixed spatial domain. Rather, this domain repre- 
sents an instantaneous statistical slice of the 3D tur- 
bulent field at the chosen location. In the present for- 
mulation, this slice is correlated (in a statistical sense) 
to the radial extent of the exhaust plume, and the 



“mean" radial scalar field at a given axial location is 
obtained by statistically averaging many realizations of 
the ID slice. The axial evolution of the exhaust plume 
is simulated by convecting this ID slice downstream 
using appropriate conservation of mass and momentum. 

Simultaneous to (but independent of) the reaction- 
diffusion simulation, turbulent mixing by the eddies in 
the flow field is implemented using a stochastic Monte 
Carlo simulation. This simulation models the effect of 
turbulent eddies on the scalar field and is captured in 
the model as a rearrangement of the scalar gradients 
using a mapping process. This mapping process (de- 
scribed in section 2 d) occurs at discrete intervals and 
can be physically correlated to the action of a turbulent 
eddy on the scalar field. The eddy size, the eddy lo- 
cation, and the mapping frequency (i.e., the turbulent 
mixing timescale) are all determined using inertial range 
scaling as applicable to full 3D turbulence. Thus, even 
though the present model is ID, the effects of the high 
Reynolds number 3D turbulence is included. It is this 
particular feature that makes the present approach ca- 
pable of capturing realistic 3D turbulence effects on the 
reaction-diffusion processes within a ID formulation. 

The ability of this model to capture realistic, high 
Reynolds number flow features in reacting shear layers, 
mixing layers, and jets has been reported in numerous 
studies in the past (e.g., Kerstein 1987, 1989, 1990; 
Menon et al. 1994). However, to study the engine ex- 
haust plumes new features had to be developed and 
incorporated. These new features are described in the 
following section. 

d. Application of LEM for jet plume studies 

Figure 1 shows a schematic of an exhaust jet plume 
and the discretization used in the ID radial LEM model. 
To incorporate the axial spreading, we assume that the 
exhaust plume is reasonably approximated by a round 
free jet. Therefore, the well-known round-jet similarity 
scaling law (Landau and Lifshitz 1959) is used to ap- 
proximate jet spreading. The resulting decrease of the 
mean jet velocity with the axial distance is given as 
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where w(0) is the axial jet velocity at the nozzle exit, d 
is the nozzle diameter, and jr 0 is the virtual origin of the 
jet (see Fig. 1). The value of the spreading rate c = 
0.11 is estimated from experimental data for a free jet. 
However, when there is a coflowing air (as in real cases 
of the engine exhaust of a flying aircraft), the jet spread- 
ing rate is decreased. This can be included by changing 
the value of c. In Eq. (1), xjd is set equal to -5 so 
that the mean flow similarity is satisfied for all x > 0. 
A constraint of the present study is that within the tur- 
bulent zone, the radial variation of the mean axial ve- 
locity m(.v, r) is neglected and. thus, u(x) is given by Eq. 

( I ) for the entire jet plume. However, the turbulent mix- 
ing process is radially dependent and is included in this 
formulation. The axial evolution of the radial line is 
related to the computational time t according to the re- 
lation dx/dt = «(*), where jr is the axial location. 


1) Molecular diffusion 

The small-scale molecular diffusion effects are ob- 
tained by numerically solving the reaction-diffusion 
equation: 


dj\ 

dt 


p 


7rP prW ' 


( 2 ) 


within the ID radial domain. Here, Y t , IV*, V*. and w k are, 
respectively, the fcth species mass fraction, molecular 
weight, diffusion velocity, and molar production rate. Also, 
p is the mixture density, and the kth species diffusion 
velocity is determined using Fick's law of diffusion: V t = 
~(D^/X t )dX t /dr, where X , and D" are, respectively, the 
*th species mole fraction and diffusion coefficients. The 
diffusion coefficients and the molar production rates for 
each species are obtained from the general-purpose CHE- 
MKIN library (e.g., Kee et al. 1989). Thus, differential 
diffusion effects are also included. 

The coordinate r in Eq. (2) is the local radial coor- 
dinate that is discretized into cells of equal size A r (see 
Fig. 1). The cell size is chosen to resolve the effects of 
all turbulent eddies present in the flow, and the domain 
is prescribed as Irl < L(x), where L(x) is the integral 
length scale defined as L(x) = c(x - x„). This relation 
for L(x) allows inclusion of the spreading rate in the 
model without actually modeling the free jet. Thus, for 
example, the increase in the jet volume (due to large- 
scale entrainment) as a function of axial distance (which 
increases the local integral length scale) is included ex- 
plicitly into the model. In contrast, in the PDF method, 
although the jet spreading rate can be enforced, the re- 
sulting mixing process cannot be physically imple- 
mented since mixing is modeled using an ad hoc model. 

The reaction-diffusion equation (2) is solved using a 
second-order, finite difference scheme (a sixth-order 


scheme was also used; however, since all length scales 
are resolved in the ID domain, it was determined that 
the computationally efficient, second-order scheme was 
sufficient) and time integrated at the local time step 
using a second-order scheme. For stability, the time step 
is determined by the minimum of the diffusion and 
chemical time steps. 

Note that in Eq. (2) the convection terms in the con- 
ventional species equations (i.e., the terms u,dYJdx,) are 
not included because turbulent transport by the velocity 
field is explicitly included as a distinct (and indepen- 
dent) mechanism, as described below. 


1 


! 


I 


2) Turbulent convection 

The effect of turbulent stirring (or convection) is rep- 
resented by stochastic mixing events that interrupt the 
deterministic solution of the reaction-diffusion equa- 
tions. Each mixing event (incorporated as a mapping 
process) is interpreted as the action of a single eddy on 
the scalar field. Three random quantities must be spec- 
ified before implementing the mapping or rearrangement 
event: 1) the segment (or eddy) size, 2) the location of 
the event, and 3) the frequency of the event. The eddy 
size, (, is determined randomly from a distribution (or 
PDF) of eddy sizes, f(C ) rj «= ( < L. The event location 
is randomly chosen within the ID domain using another 
distribution function, and the event rate (or time) is 
determined using an analogy between fluid dispersion 
and turbulent diffusivity (Kerstein 1987, 1989). Detailed 
discussions of these parameters have been reported else- 
where (Kerstein 1989; Menon et al. 1994). However, 
new modifications were needed for plume studies and 
are described below. 

The eddy size ( is chosen randomly from a power- 
law distribution f(() within the range tj < C < Ljx\ 
r), where L m {x\ r) is the largest eddy at the radial location 
r. The PDF f([ ) is determined from inertial range scal- 
ing rules (see Kerstein 1989, for detailed analysis): 


/( O = 


3Ljx; r) [(Z. m (x; r)/Tj) v ' - l]\f (x; r) 


(3) 


in terms of L„(x; r). The local eddy size L m (x\ r) is 
related to L(x) by the relation 


LJx\ r) 


L(x), 

2[L(x) - Irl], 


Irl < Z,(x)/2, 
Irl > L(x)/2. 


Equation (4) ensures that the eddies taking part in the 
mixing processes are not larger than the integral scale 
and that the chosen eddies physically reside in the jet 
shear layer. These requirements are consistent with ex- 
perimental observations of large-scale entrainment and 
small-scale turbulent mixing in jet flows. As noted ear- 
lier, this information cannot be incorporated in the PSR/ 
WMR or PDF models. 

Once the eddy size is chosen, the effect of this eddy 
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Table 1 Detailed gas phase reaction mechanisms for stratosphere (Miake-Lye cl al. 1993). The mechanism contains 24 species. 


No. 

Reactions 

A (cm' molecule 's') 

n 

E/R <K) 

1 

0 + 0 : + M=^0, + M 

O, reactions 

3 .00 x 10 » 

-2.3 

0 

■> 

O + O. => 0 , + o. 

8.00 x 10 12 

00 

2060 

3 

H + O. + M => HO. + M 

5.20 x 10 2 * 

-1.6 

0 

4 

H + O, => OH + 0, 

1.40 x 10 10 

0 

470 

5 

HO, + H =* OH + OH 

7.30 x 10 " 

0 

0 

6 

HO : + H => H O + 0 

1.60 x 10 

0 

0 

7 

HO : + H => H : + O, 

6 50 x 10 12 

0 

0 

8 

OH + 0 => H + 0, 

2 20 x 10 " 

0 

-120 

9 

HO ; + O => OH + 0. 

3 00 x 10 " 

0 

-200 

10 

H ; 0 ; + O => OH + HO. 

1.40 x 10 12 

0 

2000 

II 

OH + HO. => H O + O. 

4.80 x 10 n 

0 

-250 

12 

OH + 0, => HO. + O, 

1.60 x 10 12 

0 

940 

13 

OH + OH =» H O ♦ 6 

1.60 x 10 12 

0 

940 

14 

OH + OH + M => H O, + M 

6.60 x 10 - >y 

-0.8 

0 

15 

OH + H.O. => H O + HO. 

2.90 x 10 12 

0 

160 

16 

OH + H2 => H.O + H 

5.50 x 10 12 

0 

2000 

17 

HO. + HO. => H O. + 0. 

2.30 x 10 n 

0 

-600 

18 

HO. + 0, => OH + 20. 

1.10 x 10 14 

0 

500 

19 

N + O. => NO + O 

NO, reactions 

1.50 x 10 11 

0 

3600 

20 

N + NO =» N. + 0 

3.40 x 10 " 

0 

0 

21 

NO O + M => NO. 4- M 

4.70 x lO- 2 * 

-1.5 

0 


NO. + 0 => NO 4- O. 

6.50 x 10 12 

0 

-120 

23 

NO. + 0 4- M => NO, 4- M 

8.10 x 10 -2- 

-2 

0 

24 

NO, + 0 => NO. + 0, 

1.00 x 10 11 

0 

0 

25 

H + NO. => OH ■+■ NO 

4.00 x 10 10 

0 

340 

26 

NO 4- O, =* NO. + O, 

2.00 x 10 12 

0 

1400 

27 

NO + HO. => NO. 4- OH 

3.70 x 10 12 

0 

-250 

28 

NO + NO] => NO] + NO. 

1.50 x 10 " 

0 

-170 

29 

OH + NO + M => HONO 4- M 

1.90 x 10‘ 24 

-2.6 

0 

30 

OH + NO, + M =o HNO, 4- M 

2.20 x 10 22 

-3.2 

0 

31 

OH + HONO => H.O ^ NO, 

1.80 x 10'" 

0 

390 

32 

OH + HNO, =* H.O + NO, 

7.20 x 10 15 

0 

-785 

33 

NO. + 0, => NO, + 0, 

1.20 x I0 n 

0 

2450 

34 

NO. + NO, + M => N.O, + M 

1.00 X 10 20 

-3.9 

0 

35 

N.O, + M => NO. 4- NO, + M 

1.70 x 10 s 

-4.4 

11 080 

36 

HO. + NO, =* OH + NO. + Oj 

2.50 x 10 12 

0 

0 

37 

HO. 4- NO, => HN03 + 6, 

1.90 x 10 < 2 

0 

0 

38 

OH 4- NO, => HO, + NO, * 

2.30 x 10" 

0 

0 

39 

SO 4- 0, ^ SO. 4- o 

SO, and CO, reactions 

2.60 x 10 14 

0 

2400 

40 

so + 0, => so] 4- o. 

3.60 x 10 12 

0 

1 100 

41 

SO + NO, => SO, + NO 

1.40 x 10 " 

0 

0 

42 

HSO, + 6, => HO, 4- SO, 

1.30 X 10 12 

0 

330 

43 

SO, + OH 4- M => HSO, + M 

4.50 x lO' 2 ' 

-3.3 

0 

44 

SO 4- OH => SO, + H 

8.60 x 10 " 

0 

0 

45 

SO, 4- H.O =* H.SO, 

1.20 x lO” 

0 

0 

46 

CO + OH =* CO, + H 

1.60 x 10 n 

0 

0 


on the scalar field is implemented as a stochastic rear- 
rangement event called block inversion (Kerstein 1987) 
that mimics the effect of turbulent stirring on the scalar 
field. Given the eddy size ( from f (( ), the block inversion 
process involves first, locating randomly the center of the 
eddy within the radial line, and second, selecting the scalar 
field(s) Y k (r t t) in an interval (r 0 - (12, r 0 4- ((2) and 
replacing it with the scalar fields Y k (2r 0 - r ; r). Physically, 
inversion results in an instantaneous 1 80° rotation of the 
scalar fields contained in the eddy. Although this instan- 
taneous mapping is an artifact of the model (since turbulent 
eddies take finite time to transport the scalars), past studies 


(Kerstein 1987, 1989) have shown that this discrete map- 
ping does not cause significant errors since it occurs at the 
local mixing time, and a multitude of these mappings is 
ensemble averaged to obtain the statistically “mean” field. 
Molecular Fickian diffusion subsequently smoothens the 
sharp scalar gradients caused by these inversion events. 
In this manner, both diffusion and stirring play competitive 
and collaborative roles during the mixing process. Since 
all the length scales are resolved, no special care is needed 
during the diffusion step to deal with the sharp scalar 
gradients caused by the turbulent stirring. In an axisym- 
metric plume the cell volume must increase with radial 
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Table 2. Reduced gas phase reaction mechanism for stratosphere 
(Wang and Chen 1997). The mechanism contains 14 species. 

No. 

Reaction mechanism 

(I) 

30, = 20, 

(ID 

20H + 0. = 20, + H.O 

(III) 

20H = 0, + H, 

(IV) 

20, + NO = NO, + O, 

(V) 

OH + NO = HONO 

(VI) 

OH + NO, = HN03 

(VII) 

20, + SO. = 0, + SO, 

(VIID 

SO, + H,0 = H,SO, 

(IX) 

0, + CO = O, + CO. 


Reduced reaction rates 

(O, = 0*1 + a)3 + U) 14 - Oil 8 + a>22 + cj23 + w 25 - a)26 - a>27 

- <»)28 + a>32 - a>36 - a>37 + a>4 1 

<*>>,= w6 + all + iul3 + o)15 + o>16 4- tu31 + <d32 

a)),, = w7 — o)16 

= a)2I - a)22 - o)25 + o)26 + o)27 + o)28 + o)31 - o>4l 
= o)29 — 0)3 1 
, — o)30 — o)32 + o)37 
ok,, = o42 

"lx = o46 


distance even though a constant cell size is used to 
facilitate the mapping process. Thus, a direct imple- 
mentation of the mapping would cause mass conser- 
vation errors. Therefore, to conserve mass each mapping 
from a cell “a” with volume V a containing a scalar Y a 
to a cell “b” with volume V b containing the scalar Y bs 
the new (denoted by the superscript a) mapped scalar 
value Y n b becomes 


= [hlo - p» w t e*L}j> 

b <P- - PM + P>V b 


(5) 


Here, p a and p b are, respectively, the mass density in 
cells “a" and “b.” This mapping conserves the scalar 
fields even when the cells have different volumes. 

The second element of the turbulent stirring process 
requires locating the turbulent eddy within the local ID 
domain. The location of the block is determined ran- 
domly from a distribution function: 


1/2 

0 < lr 0 l < L(x). (6) 

Finally, the frequency of the stirring event A [or the 
timescale for the mapping, r slir = (AL) -1 ] is determined 
based on inertial range scaling laws [details are given 
in Kerstein (1987) and Menon et al. (1994) and, there- 
fore, omitted here for brevity]: 


G(lr 0 l) = 


lr 0 l 


L(x) 


V 

L(X) 

in (Ux)i - i 

ux r 

V 

1 - (v/Ux))^ 


3) Dilation process 

The x dependence of the axial component u of the 
mean velocity implies a mean flow and a radial spread- 


ing due to entrainment. These features have to be in- 
cluded to capture the global dynamics of the jet. Using 
mass conservation, the radial mean motion can be de- 
termined from the continuity equation: 


v dv _ du 

r dr dx 


( 8 ) 


Since u = u(x), Eq. (8) implies u = -( r!2)du/dx . To 
incorporate the effect of the radial motion within the 
ID formulation, the radial grid is stretched at a rate 
prescribed by Eq. (8) without changing the scalar value 
at each cell. Stretching is implemented deterministically 
by means of dilation events at regular time intervals. 
By incorporating radial flow, the algorithm conserves 
total scalar flux corresponding to the quantity ttu(x) 
j!* \r\pY k (x, r) dr in the simulation. 


e. Chemical reaction mechanisms 

The study of atmospheric chemistry within the mixing 
zone requires incorporation of many species and their 
associated chemical reactions. Depending upon the al- 
titude, ambient and engine exhaust conditions as well 
as both gas phase and heterogeneous chemical processes 
have to be included. Although detailed kinetics can be 
simulated using the LEM approach, the computational 
cost can be quite significant. Therefore, using sensitivity 
analysis (Menon and Chen 1995; Wang and Chen 1997), 
a detailed mechanism (used earlier by Miake-Lye et al. 
1993) of 24 species, 46 reactions for the lower strato- 
sphere was reduced to a 14 species, 9 reactions mech- 
anism. Tables 1 and 2 show, respectively, the detailed 
and the reduced mechanisms. Reduction of the reaction 
set is achieved by an implicit assumption that some of 
the species are under quasi-equilibrium. However, the 
accuracy of the reduced kinetics has been verified by 
using WMR and PDF calculations (Menon and Chen 
1995; Wang and Chen 1997). Since the reduced chem- 
istry is computationally very efficient, it has been em- 
ployed for the present studies. 

Heterogeneous chemical processes in the near-field 
exhaust plume can significantly modify chemical pro- 
cessing in the plume. Processes such as homogeneous 
binary nucleation, soot adsorption, condensation on aer- 
osols, and aerosol coagulation will be incorporated in 
a later study. Here, using simplified reactions occurring 
on solid condensed water (on soot particles), the effect 
of the following heterogeneous reactions on the gas 
phase species distribution (in particular, NO r , SO,, 
H 2 S0 4 ) has been studied: 

N 2 0 5 + H 2 0(c) -4 2HN0 3 (c) 
and 


S0 3 + H 2 0(c) H 2 S0 4 (c). (9) 

Here, (c) indicates the condensed (solid) state. Conden- 
sation of water on soot particles is possible because the 
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Fig. 2. (a) Variation of the scaled mean mixture fraction *£/* in the radial direction, (b) 
Variation of normalized rms fluctuation of mixture fraction CLJCLd in the radial direction 


temperature of the plume drops rapidly due to turbulent 
mixing of the exhaust plume with the ambient air (which 
is at a temperature of around 215 K). The reactions are 
assumed to occur immediately after the gas species is 
absorbed into the solid water. Thus, the heterogeneous 
reaction rates are proportional to the surface area of the 


condensed water and the flux of gas species N 2 0 5 and 
S0 3 toward the condensed water surface. Once formed, 
both nitric and sulfuric acid remain in the solid state. 
Thus, evaporation of the condensed nitric acid is not 
covered in this study. A sticking coefficient (which ac- 
counts for the fact that a collision between the gas mol- 
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Table 3. Initial conditions for the near-field plume calculations for 
the HSCT (Brown et al. 1996a). 



Exhaust flow 
condition 

Ambient condition 

Temperature (K) 

561 

205.2 

Velocity (m s _l ) 

620(450) 

0.0 

Pressure (atm) 

0.0573 

0.0573 

Mole fraction 

CO 

2.4 x 10~ 5 

20 x 10 « 

CO, 

3.2 x 10* ! 

3.5 X10 * 

H 

1.0 x 10 7 

2.8 xl0 Jl 

O 

1.6 X10-’ 

3.1 x 10 13 

h 2 

1.0 x 10'* 

5.2 x 10 7 

o 2 

0.159 

0.2098 

OH 

1.0 x 10 5 

2.8 X10 15 

ho 2 

9.9 X 10 8 

3.3 X10 13 

h 2 o 2 

2.4 X 10'* 

2.1 X 10 12 

h 2 o 

3.0 X 10 J 

4.8 X10-* 

N, 

0.779 

0.789 

NO 

4.3 X10-5 

5.9 XIO " 

NO, 

4.8 X10‘ 

4.3 x 10' 10 

HNO, 

0.0 

3.3 X10* 

SO 

0.0 

0.0 

SO ; 

6.9 X10~ 6 

0.0 

so\ 

0.0 

0.0 

HSO, 

0.0 

0.0 

h 2 so 4 

0.0 

0.0 

o 3 

0.0 

2.5 X10 * 


ecule and condensed solid water may or may not result 
in the capture of the gas molecule) of 0.1 is used to 
represent the probability of gas being absorbed in solid 
water after hitting the surface (Fukuta and Walter 1970). 
The flux of gas species to the surface is determined using 
gas kinetic theory. To model the droplet surface area 
and water condensed on active nuclei, the active nuclei 
distribution is divided into a number of groups. In each 
group, the nuclei are assumed to have the same radius. 
Based on experimental data for typical engine gas ex- 
haust (Miake-Lye et al. 1994; Brown et al. 1996a), a 
lognormal soot distribution is assumed. Water conden- 
sation is modeled by a kinetic process of aerosol droplet 
growth on preexisting nuclei particles and is based on 
a model developed by Fukuta and Walter (1970) for 
droplet growth by vapor diffusion. 

In the current implementation, all particles are as- 
sumed to travel with the mean fluid velocity, and dif- 
fusion between particles is ignored. These simplifica- 
tions are not necessary but are invoked here for com- 
putational expediency. The soot particle size distribution 
is divided into 1 1 groups. Thus, there are 1 1 more equa- 
tions to solve with the source terms given by the kinetic 
model. This is equivalent to having 1 1 more species in 
every computational cell but, unlike the gas species, 
these 1 1 groups only take part in the turbulent stirring 
and dilation processes. 

3. Numerical implementation 

The numerical implementation of LEM is relatively 
straightforward. For computational purposes, the radial 
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Table 4 Comparison of predicted mole fractions of NO and NO : along the jet centerline. 




61 m 

306 m 

613 m 

1128 

m 

* 

NO (Brown et al 1996a) 

2.87 x 10 5 

1.74 x I0-* 

8.06 x 10 7 

4.59 x 

10 7 

, 

NO (spreading rate 0 011. present) 

7.82 x 10" 

1.41 x 10 * 

7.85 x 10 7 

4.13 x 

10 7 

- 

NO (spreading rate 0.055. present) 

7.77 x 10 * 

1.77 x 10 h 

8.92 x 10 7 

4 46 x 

10 7 


NO (spreading rate O il. present) 

4.75 x 10'* 

1.48 x 10 * 

8.83 x 10 7 

4.43 x 

10 7 


NO (W'ang and Chen 1997) 

— 

3.0 x 10" 

1.6 x 10 * 





NO : (Brown et al. 1996a) 

3.47 x 10 

2.53 x 10 7 

1.24 x 10 7 

7 79 x 

10" 


NO. (spreading rate 0.01 1. present) 

8.39 x 10 7 

1.68 x 10 7 

1.00 x 10 7 

5.77 x 

10" 


NO. (spreading rate 0.055, present) 

8.71 x 10 7 

2.13 x 10 7 

1.15 X 10 7 

6.30 x 

10" 


NO. (spreading rate 0.11. present) 

5.28 x 10 7 

1.79 x 10 7 

1.13 X 10 7 

6.19 x 

10" 

1 

NO. (Wang and Chen 1997) 

— 

3.1 x 10 7 

1.8 x 10 7 

— 
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I domain -L(x) ^r< L(x) is discretized into cells small 
enough to resolve the smallest eddy in the domain. Since 
the range of eddies 17 < f < L(x) varies with the Reyn- 
olds number, the resolution must be fine enough so that 
grid effects are negligible for the entire domain of in- 
terest. For example, at x/d = 10 a resolution of 2000 
cells was sufficient to resolve all the local scales, while 
at xld = 1000 a resolution of 20 000 cells was required. 
Although this resolution requirement appears excessive, 
due to the ID nature of the simulation model, the entire 
plume with 16 species can be simulated on a workstation 
(SGI Power Indigo) in about 12 h. 

On the ID domain, the reaction-diffusion equations 
are solved deterministically at the (stability based) time 
step. The reaction-diffusion evolution is interrupted by 
the turbulent stirring process that occurs at a time in- 
terval based on the event frequency [Eq. (7)]. The mean 
flow advection is then implemented using the dilation 
process to obtain the v velocity and the stretched grid 
at the next downstream location. Reaction-diffusion and 
turbulent stirring continues at the new axial location. 
This process is continued until the downstream outflow 
■ boundary is reached. 

8 The main advantage of this method is that with proper 
p choices of the inflow and ambient conditions, the evo- 
lution of the scalar fields (with and without chemical 
& reactions) in the jet plume can be obtained as a function 
I of radial location and time, that is, Y k = Y k (r, t). Then 

I the time domain is transformed to an axial location to 

I obtain the spatial scalar field, Y k = Y k (x t r). The spatially 
t resolved statistics of the fluctuating scalar fields are ob- 
* tained by averaging a large number (greater than 50) of 


realizations (using different initial random seeds) of the 
jet plume. 

4. Results and discussion 

To validate the present formulation, simulations of 
free jets were first carried out. The effects of changing 
the Reynolds number and the Schmidt number were 
investigated and compared to experimental data. Sim- 
ulations for Re = 5000 and 20 000 and for Sc = 0.7 
(corresponding to air) and Sc = 600 (corresponding to 
liquid) were carried out using the mixture fraction ap- 
proach. Thus, the reaction-diffusion equation [Eq. (3)] 
was replaced by a single equation for the mixture frac- 
tion: d£ldt = Dd 2 £/dr 2 . Here, £(x, r) is the mixture frac- 
tion and D is the diffusion coefficient. 

Since the jet is evolving in the self-similar regime, 
the mixture fraction field should also show self-simi- 
larity. This is demonstrated in Fig. 2a, which shows the 
mean mixture fraction radial profiles in the similarity 
form g(x, r) = kg{r))f\ at various axial locations plotted 
as a function of the scaled radial coordinate 77 = r!(x 
— x 0 ). Here, g(r j) is the similarity form of the mixture 
fraction and \ “ (* “ x 0 )/d. The figure shows that self- 
similarity is achieved in the present calculation and that 
the result is consistent with experimental data (Dowling 
1988). The radial profile of root-mean-square (rms) fluc- 
tuation of mixture fraction, shown in Fig. 2b, also ex- 
hibits self-similarity as expected. The occurrence of the 
peak rms fluctuation, where the slope of the mean profile 
is greatest, is a general feature in free jets and has been 
reproduced by the present model. The model predictions 


f 

f 

f 


Table 5. Comparison of predicted temperature and species mole fractions at a fixed axial location. 


306 m 

Brown et al. 
(1996a) 

Present study 
with 

molecular diffusion 

Present study 
without 

molecular diffusion 

Present study 
with molecular diffusion 
and heterogeneous reaction 

Temp 

229.7 K 

213.9 K 

213.9 K 

217.2 K 

OH 

1.54 X 10 9 

1.06 X 10 10 

0.87 x 10'° 

1.29 X I0 >° 

NO 

1.74 x 10" 

1.77 X 10" 

1.04 X 10" 

1.39 x 10" 

NO, 

2.53 X 10 7 

2.13 x 10 7 

1.24 x 10 7 

1.67 X 10 7 

HNO, 

3.98 X 10" 

2.80 X 10" 

2.62 X 10" 

3.70 X 10" 

S0 2 

3.05 X JO" 

2.37 x 10- 7 

2.30 X 10 7 

3.18 x 10" 

SO3 

1.02 X 10" 

1.01 X 10 10 

0.96 x 10 '° 

0.29 x 10 10 

h,so 4 

3.19 X 10" 

5.76 x 10" 

5.70 x 10" 

1.22 X 10" 
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Fig. 4. (a) Axial NO, response variation with spreading rate for an initial jet exit speed of 620 
ms' 1 , (b) Axial O* response variation with spreading rate for an initial jet exit speed of 620 m 
s"L (c) The 0 } depletion along the exhaust centerline (radially integrated). 


deviate from the data near the edges of the jet. This is 
due to the intermittency (in the experiments) that was 
not included in the present formulation. 

Three different realizations of the scalar field are 
shown in Fig. 3 for Re « 5000, Sc « 0.7, and at a 
location of x - (■* - x Q )fd = 25. This is quite consistent 
with the experimental observations of Dahm and Di- 
motakis (1987). The scalar mixture fraction fluctuates 


about a constant scalar level near the centerline with 
a relatively sharp drop at the boundaries of the shear 
layer. This behavior is seen both in water shear layers 
(Sc = 600) and in air shear layers (Cruyningen et al. 
1989). 

The LEM approach was then implemented to study 
aircraft jet plumes. Simulations were carried out for rep* 
resentative conditions for the HSCT and are summarized 
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(C) 



in Table 3. To include the effects of the ambient flow 
surrounding the jet (which is not included in the present 
formulation but would effectively decrease the spreading 
rate), the spreading rate parameter c and the jet exit ve- 
locity were varied. The domain simulated was a region 
from the engine exit to 1000 jet diameters downstream. 
For a reference initial jet diameter of 1.2 m, this implies 
a domain of 1 .2 km. Although this is an appropriate region 
for the near-field plume study (Quackenbush et al. 1993; 
Brown et al. 1996a), it should be noted that by the end 
of this domain, the engine exhaust plume could be en- 
trained into the wingtip vortex. The plume-vortex inter- 
action is not included in this study. However, the present 
simulation provides an appropriate boundary condition for 
the full 3D simulation of the plume-vortex interaction 
process (Menon and Wu 1997). 

Several quantities describing the depletion of NO, and 
0 3 , such as the responses of NO, N0 2 , 0 3 , and NO x 
(sum of NO and N0 2 ) and the depletion of 0 3 with 
respect to the initial NO x , have been computed for all 
the simulations. However, only representative results are 
shown here. The species response is defined as 


Tine with 
the shear 
it layers * 
__en et al. ' 

to study 
for rcp- 

-nimarizcd 


R, = 


m,(t) ~ m, , 


m, 


( 10 ) 


where m,(r) is the mass flux of species i evaluated over 
the cross section of the jet at a certain time that then 
corresponds to a certain downstream location of the jet. 
Also, m, lotll is the sum of the mass flux of the ith species 
originating from the nozzle exit, and the mass 

flux entrained into the jet, m, enlnMMd . A response R, = 0 
implies no change of the total amount of the species. 
Negative or positive responses correspond to destruction 


t; 

K 


* 


or creation of the zth species, respectively. The depletion 
of 0 } is defined relative to the amount of NO, origi- 
nating from the nozzle exit, m NO , and is given as 
[O3Lepiet.cn = [m 0j (r) - wt 0enuained ]/m NO( . By definition, 
this quantity is always negative and gives an estimate 
of the amount of 0 3 destroyed by NO, in the jet plume. 

Species concentrations along the axis of the jet plume 
at different locations are compared to the results from 
Brown et al. (1996a) and Wang and Chen (1997) in 
Tables 4 and 5. The results for NO and NO : show rea- 
sonable agreement after x/d > 250 (greater than 300 m) 
with earlier calculations. The differences in the near 
field, that is, x/d < 50 (61 m), may be due to the non- 
similar initial evolution of the jet immediately down- 
stream of the engine exit plane. This feature is not cap- 
tured since the present scaling requires jet similarity. 
The spreading rate does not have much effect on the 
axial species concentrations of NO and N0 2 because 
the initial NO and N0 2 are the same and no ambient 
NO and N0 2 are entrained. Although the spreading rate 
does show some effect on the data in the very near field, 
by x/d > 500 (greater than 600 m), its effect has di- 
minished considerably. 

Some results for other radicals are also shown in Table 
5. Only limited comparison was possible with other 
model predictions and, therefore, no major conclusions 
can be made at this time. However, it is evident that 
inclusion of molecular diffusion and heterogeneous ki- 
netics does have an impact on the OH radical and sul- 
furic acid concentrations. The changes in OH concen- 
tration have an impact on ozone depletion, as discussed 
below. The decrease in the gas phase H 2 S0 4 formation 
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Fig. 5. (a) Comparison of axial variation of NO, response predicted by the present LEM 
simulation and PDF simulation (Wang and Chen 1997). (b) Comparison of axial variation of O. 
depletion predicted by the present LEM simulation and PDF simulation (Wang and Chen 1997). 


is due to the heterogeneous conversion of SO, to the 
condensed form of H,SO„ on soot particles. 

Figure 4a shows the NO, response in the jet plume. 
Initially, the NO, response has a sudden drop (a 5%- 
6% decrease) and then reaches an almost constant value. 
Notice that the response does not increase consistently 
with decrease in spreading rate. This can be shown to 
be related to the competition between various reaction 


mechanisms as the ozone concentration changes (de- 
creases with decrease in spreading rate). The reasons 
for the initial sudden drop are the reaction (Rl), NO + 
0 } = NOj + O,, due to entrainment of O, from ambient 
air, and the reaction of converting NO, to nitric acid 
(R2), NO, + OH + M = HN0 3 + M. Initially, the 
concentration of NO is much larger than N0 2 , and the 
destruction rate of NO depends on the concentration of 
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0 3 entrained from the ambient air (R 1 ). Since the ozone 
entrainment increases with spreading rate, it can cause 
an effective decrease in OH (via reactions 11, 12, and 
18; see Table 1). This will decrease the conversion of 
NO : to nitric acid (R2), thereby decreasing the destruc- 
tion of NO,. However, if the ozone concentration in- 
creases above a certain level, Rl dominates, leading to 
a large production of NO : that is then converted to nitric 
acid by R2. The competition between these reaction 
pathways leads to the NO, variation with spreading rate. 
It is seen in all cases that, after around x/d = 100, the 
depletion of NO is balanced by the creation of N0 2 . 

Figure 4b shows the corresponding 0 3 response. Ini- 
tially, the 0 3 response also drops from a positive to a 
negative value and then tends to almost zero. The spread- 
ing rate has an effect on the 0 3 response with the largest 
O, response corresponding to the smallest spreading rate. 
This is reasonable since the spreading rate controls the 
entrainment of 0 3 . The initial drop is due to the initial 
rapid entrainment of 0 3 , which allows the reaction NO 
+ O, = NO : + 0 ; to deplete O v The initial variation 
in the 0 3 response is due to the differences in the 0 3 
entrained (which increases with spreading rate) and the 
depletion of 0 3 (which depends on the NO, from the 
engine and the temperature). Thus, for the lowest spread- 
ing rate (smallest amount of 0 3 entrained), the reactions 
between O, and NO, are very effective in depleting 0 3 
since very little mixing occurs (which results in high 
temperature and high NO, concentration). However, 
when the spreading rate is increased and more ozone is 
entrained, the increased mixing in the jet enhances the 
dilution of NO, concentration and the temperature also 
decreases. As a result, the ozone response is very small. 
Further downstream, although 0 3 entrainment linearly 
increases with the distance, the radicals from the engine 
have been almost exhausted and results in almost neg- 
ligible depletion of 0 3 . 

Figure 4c shows the ozone depletion as a function of 
spreading rate and initial jet velocity. In all cases, ozone 
is depleted by around 2%-3% in the jet plume. Increase 
in the spreading rate (for a fixed jet velocity) slightly 
increases ozone depletion. There is, however, a much 
stronger dependence on the jet velocity with more ozone 
depleted when the jet velocity is reduced. This is un- 
derstandable since reducing the jet velocity implies re- 
duced convective motion (increased large-eddy turnover 
time). This allows more time for the local reaction- 
diffusion and the turbulent mixing processes to complete 
before downstream convection occurs. 

A comparison between the present results and the 
PDF jet simulation by Wang and Chen (1997) was car- 
ried out for a fixed spreading rate of c = 0.055 and an 
initial jet velocity of 620 m s _I . Figures 5a and 5b show, 
respectively, the NO, response and ozone depletion pre- 
dicted by the two methods. Although the same trend is 
seen, this study shows a reduced NO, response and a 
reduced ozone depletion in the jet plume. Around 15% 
less NO, and 30% less 0 3 are depleted when microscale 


turbulent mixing and molecular diffusion effects are in- 
cluded. Earlier, comparisons between WMR and PDF 
predictions were carried out (Wang and Chen 1997). 
Results showed that a similar decrease in depletion oc- 
curs even when mixing was modeled (using an ad hoc 
mixing timescale model). However, the overall response 
was not very high. The WMR and PDF methods pre- 
dicted, respectively, a 6% and a 5% decrease in NO, in 
the plume while the present results suggest an even 
smaller decrease: approximately 4%. Thus, it appears 
that taking into account the local turbulent mixing and 
diffusion timescales results in a reduction in NO, de- 
pletion. A similar result for ozone depletion is shown 
in Fig. 5b and indicates that the inclusion of local mixing 
effects result in an overall 0 3 response of around 0.25%, 
compared to around 0.35% as predicted by the PDF 
method. 

To determine the importance of molecular diffusion (and 
differential diffusion), a simulation was carried out by 
suppressing molecular diffusion but allowing the small- 
scale turbulent stirring to occur unchanged. This result, 
also shown in Fig. 5, suggests that accounting for turbulent 
mixing at all the spatial and temporal small scales has a 
much larger impact than including molecular diffusion. 

When heterogeneous reactions are included, NO, is 
transformed to its relatively stable form in the reservoir 
specie nitric acid. This increases the destruction rate of 
NO, (by around 15%), as shown in Fig. 6a (some results 
for other radicals are shown in Table 5). Note that most 
of the NO, destruction occurs within the first 100 jet 
diameters, and after that remains almost constant. De- 
struction of O, (not shown) is also decreased, although 
at a less significant level than NO, destruction because 
there is less NO,. However, since 0 ; is also destroyed 
by radicals such as OH and HO : , the final 0 3 response 
is almost unchanged when heterogeneous reactions are 
included. This result further emphasizes the importance 
of properly estimating the concentration of other radi- 
cals such as OH and H0 2 . 

The competition between the gas phase H^S0 4 and 
the condensed phase H : S0 4 formation reactions for the 
available S0 3 was also studied. Figure 6b shows the 
centerline SO, concentration for these cases. The S0 3 
concentration is significantly decreased when the gas 
phase H 2 S0 4 formation reaction is included. Adding the 
heterogeneous process of H 2 S0 4 condensation on soot 
particles causes a further decrease in the SO ? concen- 
tration. This behavior is reflected in the gas phase H 2 S0 4 
centerline mass fractions for these two cases, as shown 
in Fig. 6c. When condensed H 2 S0 4 is formed, the gas 
phase H 2 S0 4 concentration is significantly decreased 
since a substantial amount of S0 3 is lost from the gas 
phase to the condensed phase reaction. Note that this 
study did not address the other heterogeneous processes 
related to sulfur particles in the plume. Thus, these re- 
sults are likely to be modified when binary homogenous 
nucleation and other aerosol-related kinetics are includ- 
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Fig. 6. (a) Effect of heterogeneous kinetics on the NO, response in the exhaust plume. (b) 
Centerline variation of SO, with and without sulfuric acid formation, (c) Centerline variation 
of H,S0 4 . 


ed. Extension of this model to factor in these hetero- 
geneous processes is currently under way. 

5. Conclusions 

In this study, fluid mixing and chemical reactions in 
the near field of an engine exhaust jet plume were in- 
vestigated using a mixing mode! that explicitly incor- 


porates the large-scale mixing process as entrainment 
events and stochastically simulates small-scale turbulent 
mixing and molecular diffusion effects. A reduced re- 
action mechanism that accurately reproduces results of 
a detailed mechanism for the lower- stratosphere chem- t 

istry was used to simulate the exhaust plume dynamics j 

of a typical HSCT aircraft. i 

This study shows that, compared to predictions based L 
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on models that do not include small-scale mixing and/ 
or molecular diffusion effects, 30% less 0 3 and 15% 
less NO, are depleted in the near field of the plume. 
This suggests that the effects of the local lack of mixing 
have an inhibiting effect on ozone depletion in the near- 
field plume. Inclusion of heterogeneous kinetics in- 
volving formation of condensed nitric and sulfuric acid 
on water condensed on soot particles showed that 15% 
of the available NO, is converted into its inactive form 
but has negligible effects on O, concentration. The pre- 
dictions are also shown to be dependent on the avail- 
ability of other radicals such as OH. Although these 
results are in reasonable agreement with other predic- 
tions using less sophisticated (or no) mixing models, 
the importance of the small-scale mixing and diffusion 
processes on the near-field plume dynamics can only be 
addressed using a physically consistent simulation mod- 
el. The development of this model offers such a capa- 
bility. The impact of small-scale mixing on homoge- 
neous nucleation of sulfuric acid aerosol particles has 
not yet been addressed but is currently being incorpo- 
rated and will be reported soon. Another possible role 
of this model is to provide appropriate conditions for 
the plume-wingtip vortex interaction phase that im- 
mediately follows the plume stage. More complex 
chemistry, especially heterogeneous kinetics that occur 
in this stage, depends critically on the chemical state of 
the plume during the interaction phase. Accurate esti- 
mate of the species concentrations and aerosol particle 
distribution at the end of the plume stage will allow a 
more accurate evaluation of the plume-vortex dynam- 
os. To evaluate this capability, simulations are currently 
under way (Menon and Wu 1997) of the plume-vortex 


interaction phase using data from this study. These re- 
sults will be reported in the future. 
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in near-field aircraft plumes 
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Abstract. The potential impact of aircraft exhaust gas on entrained ambient air in the strato- 
sphere has been studied by using several numerical models to manifest the importance of mix- 
ing processes occurring at different scales. These models include perfectly stirred reactor 
(PSR), transient well-mixed reactor (WMR), partially stirred reactor (PaSR), and joint scalar 
probability density function (PDF) simulation for the near-field jet exhaust. Sensitivities of 
major species to chemical reactions are presented to illustrate the chemical transformation in 
the near-field jet regime. An accurate reduced chemical mechanism that captures major reac- 
tion pathways has been systematically developed and tested. The predicted results with the re- 
duced mechanism are found in good agreement with those from the detailed mechanism, but the 
computing time is decreased by a factor of 3. The effect of microscale turbulent mixing and 
chemistry interaction on species evolution is studied, and its impact is found most noticeable in 
the near field where the timescale of chemistry for highly reactive radicals, such as OH and 
HO:, is comparable or faster compared with that of microscale mixing. When microscale 
mixing is included, the predicted H 2 S0 4 level increases at locations where sulfate aerosol for- 
mation may be important in the near field. 


Introduction 

Engine emissions from the rapidly growing subsonic air- 
crafts and the planned high-flying supersonic aircrafts may 
significantly alter concentrations of ozone, water vapor, and 
aerosol in upper atmosphere and may lead to climate changes. 
In particular, potential ozone depletion caused by aircraft 
emissions in the stratosphere is of great concern. Conse- 
quently, the impact of emissions from the planned supersonic 
commercial aircrafts on atmospheric chemistry needs to be as- 
sessed as thoroughly as possible. Currently, several major re- 
search efforts are undertaken within NASA’s High-Speed Re- 
search Program (HSRP) to address various issues related to 
the effects of a fleet of high-speed civil transport (HSCT) 
flights on the stratospheric chemistry. Research tasks on 
similar subjects for the growing subsomc aircrafts are also 
being conducted within NASA’s Advanced Subsonic Tech- 
nology Program (ASTP). A review of these studies [Stolarski 
and Wesoky , 1995] clearly demonstrates the complexity of the 
problems and the associated uncertainties. Numerical calcula- 
tions on near-field jet plume [Kiiake-Lye et ai, 1993a, b; Kdr - 
cher, 1994, 1996, Anderson et ai , 1996; Brown et ai , 1 996] 
and experiments in the wake [Arnold et ai, 1992; Fahey et 
ai, 1995] have been implemented to quantify levels of major 
gas phase species and exhaust aerosols. Although the 
near- field jet regime lasts only a few seconds, major conver- 
sion of gas phase NO, and SO, into HN0 3 and H 2 S0 4 driven 
by hydroxyl radical occurs in the near field. Nucleation of 
sulfate aerosol particles and their possible interaction with ac- 
tive soot particles are believed to occur primarily in this re- 
gime. Modeling of near-field jet is considered an important 
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task to provide proper initial conditions of NO, and aerosol 
contents for the subsequent study on chemistry' and micro- 
physics in the intermediate vortex interaction and far-away 
dispersion regimes. One area that lacks fundamental under- 
standing is the relative importance of microscale turbulence 
and chemistry interaction in the near-field aircraft plume and 
its influence on chemical kinetics when combined with large- 
scale entrainment process. 

Since the fluid mechanics becomes complicated m the re- 
gion where the exhaust jet plume and wing vortex begin to 
interact, large-eddy simulations may be necessary for explor- 
ing the role of wing vortex on mixing process The conse- 
quential impact of wing vortex on chemical kinetics can be 
studied provided that simplified chemistry is available for de- 
scribing the key reaction paths. In addition, aerosol dynamics, 
e g., binary nucleation of sulfuric acid and water vapor, con- 
densation, coagulation, are considered important processes m 
the near-field plume. Inclusion of these processes, while tak- 
ing into account the interaction between microscale turbulence 
and chemistry, would demand substantial computing time. 
The main objectives of this study are to identify and charac- 
terize the distribution of key species in the near-field plume 
regime, to develop an accurate reduced chemical mechanism 
for gas phase chemical reactions, and to explore the impact of 
microscale mixing on the chemistry in the exhaust plume. 

Numerical Models 

Table 1 lists the major features of various numerical mod- 
els that are used to assess the relative importance of mixing 
processes at different scales on the evolution of plume chem- 
istry in the near field of exhaust plume. The macroscale mix- 
ing process refers to turbulent convection and stirring across 
the jet plume. The dominant role of macroscale mixing is to 
bring fluids from different parts of flow over a distance much 
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Table 1. Features of Various Numerical Models for Modeling Jet Plume Chemistry 


Model Assumptions PSR WMR PaSR PDF Jet With Each Cel! PDF Jet With Modified Curl’s Mixing 

Well Mixed Model For Each Cell 


Microscale well mixed 
Macroscale well mixed 
Transient effect 


PSR, perfectly stirred reactor, WMR, well-mixed reactor, PaSR, partially stirred reactor, PDF, probability density function^ Micro- 
scale well mixed means micromixing time is much smaller than the chemistry time, so microscale turbulence effect is ignored. Mac- 
roscale well mixed means macroscale turbulent diffusion time is much smaller than the residence time, so properties are spatially uni- 
form (homogeneous) It is assumed that in PSR and PaSR, there is a stirrer which makes the reactor homogeneous and either well 
mixed or not well mixed in microscale. WMR is similar to a Lagrangian box model, and PDF jet with each cell well mixed is equiva- 
lent to the standard plume flow field code [ Dash et al., 1990] and K&rchet*s [1996] 2-D jet model. 


larger than the scale of molecular diffusion and to create sub- 
stantially large interfacial areas among fluids. The microscale 
mixing refers to the turbulent diffusion process occurring at 
the microscale which is strongly enhanced by macroscale 
mixing. Only after fluids are well mixed at the microscale will 
chemical reactions take place. The resolution of numerical 
simulation for the mixing process is limited by the grid size 
which is usually much larger than the scale of molecular dif- 
fusion. Consequently, modeling of turbulence and mixing pro- 
cess for fluids within grid cells is necessary. A microscale 
well-mixed state can be assumed when the turbulent mixing 
within a grid cell is faster than chemical kinetics. 

If the mixing processes at both macroscale and microscale 
are assumed infinitely fast, the jet plume can be modeled by 
the well-mixed reactor (WMR) which is similar to the La- 
grangian box model [e g., Kiiake-Lye et al 1993b; Karol et 
al. y 1994; Karcher, 1995]. A perfectly stirred reactor (PSR) is 
a simplified version of WMR as it computes only steady state 
solutions of chemical kinetics for a mixture of different fluids 
brought together into the reactor for a period of residence time. 
PSR calculations facilitate detailed sensitivity analysis which 
is useful for identifying the major chemical reaction pathways. 

When the fast mixing assumption at the microscale is re- 
laxed, a partially stirred reactor (PaSR) is developed to in- 
clude the effect of microscale mixing. The unmixed nature is 
described by the joint scalar probability density function 
(PDF) for reactive species. To account for the interaction of 
fluid dynamics and chemistry that causes non-uniform fluid 
compositions inside the jet plume, two-dimensional modeling 
of jet is necessary to resolve the spatial distribution. For in- 
stance, the standard plume flow (SPF-II) field model devel- 
oped by Dash et al. [1990] can be used to model jet develop- 
ment with the assumption of fast microscale mixing within 
grid cells. To explore the potential impact on chemical kinet- 
ics due to microscale mixing, the joint scalar PDF is used to 
describe the unmixed nature of fluids within each grid cell. 
The PDF jet model will be computed with/without microscale 
mixing models to provide information on the relative impor- 
tance of microscale turbulent mixing and its interaction with 
chemical kinetics. Described below are details of these nu- 
merical models. 

Perfectly Stirred Reactor (PSR) Model 

The perfectly stirred reactor (PSR) is a homogeneous reac- 
tor which has a continuous mass flow in/out of the reactor as 
shown in Figure 1 [Glarborg et al., 1986]. When the fluids 


enter the reactor, it is assumed that the fluids are perfectly J 

mixed with the existing fluids in the reactor immediately. This } 

assumption is justified when the characteristic time of turbu- j 

lent mixing in both macroscale and microscale is much 
smaller than that of chemistry, so the rate of chemical process | 

is controlled solely by chemical kinetics. Turbulence effects 
are not accounted for in this model because of the assumption 
of infinitely fast mixing. Detailed chemical kinetics can be 
studied by using mean thermochemical properties. 

The conservation equation for reactive species at the steady 
state is given by 

m(Y k -Y' k )-c* k W k V = 0, (1) 

and the corresponding energy equation is described by 

mi(Y k h k -r* k h' k ) + Q=0. ( 2 ) 

A=1 

In these two equations, Y k is the mass fraction of the k\h spe- 
cies (there are K species); W k is the molecular weight of the 
fah species; V is the reactor volume; m is the mass flow rate 
through the reactor; Q) k is the molar production rate of Mi 
species per unit volume by chemical reactions; /i* is the spe- 
cific enthalpy of the Mi species; and Q is the reactor heat loss. j 

The superscript asterisk indicates the inlet conditions. The j 

characteristic flow time is residence time r = pV tm A hy- j 

brid Newton/time-integration algorithm is used to solve for 
the solution, and details of the numerical method are de- f 

scribed by Glarborg et al . [1986]. PSR is used in the current 



Figure 1 . Schematic of perfectly stirred reactor (PSR) show- 
ing the characteristics of the reactor, where rh is the mass 
flow rate, k denotes species k, and superscript asterisk indi- 
cates the inlet condition. The key parameter is the residence 
time r = pV I rh . 
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study primarily for developing reduced chemistry and for 
identifying controlling steps in 0 3 and NO, chemistry through 
sensitivity analysis. To assess different treatments for the 
mixing processes occurring in the near field of the jet plume, 
we have used more sophisticated numerical models described 
below. 


Transient Well-Mixed Reactor (WMR) Model 


The transient WMR model facilitates study of the time 
evolution of chemical kinetics for engine exhaust gas subject 
to injection of ambient air. A schematic diagram for the WMR 
is sketched in Figure 2. Similar to the PSR, a global residence 
time can be defined as 



(3) 


where A/* is the mass of gas inside the reactor, and m a is mass 
entrainment rate. Using z> as a parameter, conservation 
equations of mass, species, and energy can be written as 


dt 


_1_ 

*R 


dM R ki R 

(4) 

dt t r 

jm 1,. 
PR 

(5) 


^-7^— [Z(V/ -h RJ )Y aJ ]- 

dt CpR^R J = t J=\ PrCpr 


CprMr 

( 6 ) 


where the subscript R represents quantities inside the reactor 
and those entrained are denoted by the subscript a The WMR 
is applied to study the chemical kinetics of the engine exhaust 
gas that is continuously mixed with ambient air. For constant 
temperature cases the WMR model is similar to the Lagran- 
gian box model [Miake-Lye et al. y 1993b; Karol et ai , 1994; 
Kdrcher , 1995] when the following relation is used to bridge 
the two models: 


m a = K(t) 
Mr A * 


( 7 ) 


where A is the cross-section area of a plume and K(t) is the 
growth rate of the plume cross-section area m a Lagrangian 
box model equation 



Figure 2. Schematic diagram of a transient well-stirred reac- 
tor (WMR) showing ambient air entrained into the reactor 
with a rate of m a , where M R is the mass of gas in the reactor, 
and it increases with time. 


^ = V - (8) 

at A 

where n h n? are the species concentrations inside the box and 
the air, respectively; P x and L, are the production and destruc- 
tion rate by chemical reactions; R t is the source or sink due to 
other mechanisms. In the WMR model reported here, the air 
injection rate is determined from the mass entrainment rate 
into the jet estimated by the PDF jet plume model to be dis- 
cussed in the later session. 


Partially Stirred Reactor (PaSR) Model 

When the characteristic time of turbulent microscale mix- 
ing is not fast compared to that of chemical reaction, the im- 
pact of microscale mixing may become important. In the cur- 
rent PaSR model, the mean properties are assumed to be spa- 
tially homogeneous (well mixed in macroscale turbulence) but 
not completely mixed at the molecular level. Compared to 
WMR, the effect of microscale mixing is included in PaSR, so 
the scalar distribution inside a reactor needs to be modeled to 
manifest the micromixing process. The joint scalar probability 
density function (PDF) is used to describe the distribution of 
scalars, and therefore it contains information about the mixing 
state. The time evolution equation for the joint PDF for the 
PaSR has been derived by Chen [ 1 996] as 


- ffi 0 = - £ 

d i a =[dyf a l 


1 




k 

- z 


1 dWadVp 




(9) 


where the mean residence time is defined as r R - M R l m a , 
scalar ^ represents the set of species mass fractions and en- 
thalpy, is the set of independent variables corresponding to 
the dependent scalar variables P+iysj) is the joint PDF of 
a tilde denotes the density-weighted average, €<# is the 
scalar dissipation rate given by €<# = where T is 

molecular diffusivity, <e a# \ # = y/> represents the average 
scalar dissipation rate conditioned on a certain set of scalars, 
and S a is the chemical source term for scalar ^ The first two 
terms on the right-hand side represent the effects of chemical 
reaction and the throughflow on the joint scalar PDF, respec- 
tively. Both of these terms can be calculated directly and do 
not need modeling. The last term stands for the effect of mi- 
croscale mixing on the PDF which requires closure approxi- 
mations. PaSR includes chemical kinetics and the additional 
interactions with microscale turbulence under the assumption 
of homogenous flow (no macroscale turbulence) and further 
ignoring velocity fluctuations. The simplicity of this model 
permits us to cany out simulations with detailed chemistry 
without significant computing burdens. More importantly, 
PaSR can be viewed as a single representative grid cell em- 
bedded in a large computation scheme. Therefore PaSR is an 
ideal test bed for exploring the influence of the microscale 
mixing on chemical kinetics, for evaluating the performances 
of reduced chemistry, and for assessing current and emerging 
mixing models. 
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In this study, the effect of microscale mixing on the PDF is 
modeled by the modified Curl’s mixing model [Janika et 
al. 1979], 

- I - - — [< £afl\t = V'> P*(¥J)] 

1 4 ff?*(V'V)^(V'V)tf (¥''.¥'* I v)diifd\y- P^.O j. 

( 10 ) 

where H is the transitional probability defined as 
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marching downstream algorithm, and the details can be found 
in the works of Chen et al. [1987, 1989]. The modified Curl's 
mixing model (equation (10)) was used to simulate the effect 
of microscale mixing process. 

The PDF numerical model was validated by comparing 
computed results without the effect of microscale mixing with 
those obtained from the SPF-II code by Brown et al. [1996] 
using the same mechanism and initial conditions for a HSCT 
case. The computed centerline mean concentrations of NO x , 
NO* SO 2 , and H 2 S0 4 are found in good agreement despite 
many differences in the numerical models In addition, both 
models predict the NO jr /NO > ratio in excess of 80% in accord 
with the measurements by Fahey et al [1995] from the 
NASA ER-2 high-altitude aircraft. Furthermore, we com- 


O 1 ) 

The amount of mixing is modeled through a characteristic 
mixing time which is related to the turbulence timescale for 
micromixing, t mLX ~k Is , where k is turbulence kinetic en- 
ergy and e is turbulence dissipation rate. Because of the 
multidimensional nature of equation (9), Monte Carlo tech- 
nique is used to solve the PDF evolution equation. In the con- 
text of simulating jet plume using a Pa SR, the air injection 
rate, similar to that in the WMR, is needed, and its value is 
estimated from the PDF jet plume model described next. 

Probability Density Function (PDF) Simulation for 
Near-Field Jet Exhaust 

A proper treatment of the interaction between turbulence 
and chemical kinetics is important to the predictions of ozone 
and NO, changes in the engine exhaust gas. The standard 
plume flow (SPF-II) field model by Dash et al. [1990] as- 
sumes that the gas inside each computational cell is perfectly 
mixed. The microscale mixing effect is therefore neglected in 
the SPF-II code, and its impact on chemical kinetics has not 
been assessed. To account for the microscale mixing, the joint 
scalar PDF is solved for by the following transport equation 
[Pope, 1985] using the Monte Carlo technique, 


* 

E 


* 

E 


^ [p < e a0 I * = v > 

0=1. /j=i dy a oy 

( 12 ) 

where the turbulent transport term can be modeled by a gradi- 
ent diffusion model using the turbulence time and fluctuating 
velocities from the Reynolds stress model, 

<u’ a \t=<r>p i =-c s - »y, (B) 

e acp 

where C s is an adjustable constant (= 0.25 in jet flows), 
uj* = u k - u k is the density -weighted velocity fluctuation, and 
<u k "\4> = y> is the average fluctuation conditioned on a cer- 
tain set of scalars. The flow field is solved by a Reynolds 
stress turbulence closure model, which directly solves the 
modeled transport equations for Reynolds stress, with a 



x(m) 


Figure 3. Centerline mole fractions of major nitrogen, sulfur, 
and hydrogen species for B747 jet engine. They compare well 
with Karcher's [1996] results (not shown). 
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pared our predicted results with Kdrcher's [1996] B747 case 
using the same chemical mechanism and initial conditions. 
Figure 3 presents our predicted centerline mole fractions for 
NO,, HO„ and SO,, which agree with Kircher's numerical re- 
sults (not shown). 

Chemical Reaction Mechanisms 

Detailed Gas Phase Mechanism 

Tables 2a and 2b list a detailed gas phase mechanism con- 
taining 24 species with 60 elementary reaction steps. It in- 
cludes reactions of species O,, HO„ NO,, SO,, and CO, that 
are important in the near-filed jet plume in the stratosphere. 
The reaction rates are basically taken from data presented by 
Karcher [1996], which cover reaction rates below 1000 K and 
most of which represent linear or nonlinear fits to a large 
number of the data by National Institute of Standards and 


Table 2a. Gas Phase Reaction Mechanism: Binary Reactions 


Reaction Rate Constant 

(cm 3 molecule’ 1 s’ 1 ) 


(Rl) 0 + 0,-+ 20 , 

(R2) H + Oj -> OH + O, 

(R3) H + OH-*0 + H 2 
(R4) H + HO? -+ 20H 
(R5) H +■ HO, -► H, + O, 

(R6) H + HOj -> O + H 2 O 
(R7) OH + O -► H + O, 

(R8) OH + O, -► HCb + O, 

(R9) OH + Hj-»H+ H,0 
(R10) OH + OH ->Ot H,0 
(Rl 1) OH + HCb -* Oj + HjO 
(R12) OH + H:0: -> HO: + H:0 
(R13) HO 2 + O -> OH + O: 
(R14) HO 2 O, -» OH + 20: 
(R15) HO 2 + HO: -» H,0, + O 2 
(R16) H,0, + O -» OH + HO, 
(R17) H,0, 1 - H -» OH + H,0 
(R18) H ,0, + H -► HO, + H: 

(R19) NO + O, -♦ NO, + O, 
(R20) NO + HO: -* OH + NO 2 
(R21) NO + NO) -» 2 NO 2 
(R22) NO: + O -» NO + O, 

(R23) NO: + O, -» NO, + O: 
(R24) NO, + H -» OH + NO 
(R25) NO, + NO, -* NO + NO, + 
(R26) NO, + O -» NO, + O, 
(R27) NO, + OH -► HO, + NO, 
(R28) HNO, + O -► OH + NO, 
(R29) HNO, + H -> NO, + H, 
(R30) HNO, + OH -» NO, + H,0 
(R31) HNO, + O-* OH + NO, 
(R32) HNO, + OH -» NO, + H 2 O 

(R33) SO + O, -> O + SO, 

(R34) SO + O, -► SO, + O, 

(R35) SO + OH -► H + SO, 

(R36) SO + NO, -+ NO + SO, 
(R37) SO, + O, -> O, + SO, 
(R38) SO, + O -► O 2 + SO, 

(R39) SO, + H,0 -» H,SO, 

(R40) HSO, + O, -► HO, + SO, 

(R41 ) CO + OH -♦ H + CO, 


1.21 x 10’"exp(-2I25/T) 

1.15 x 10 lo exp(-436/T) 

8 10 x j o’ 21 t 2 ’ exp( 1 950/T) 

2.8 x 10’ 10 exp(-440/T) 

6.9 x 10’"exp(-636.97r) 

2.80 x l0' 12 T o44 exp(-677.9/T) 
1.83 x 10’" exp(173.3/T) 

1.9 x I0' l2 exp(-1000/T) 

lit X 10’ l4 T l ” exp(-15897T) 
8.34 x 10’ 17 T‘ 54 exp(355.7/T) 

5.09 x I0'"exp(72.6/T) 

2.13 x 10’ 13 T u exp(-179.8/T) 
2.71 x 10’"exp(-224/T) 

14 x 10' l4 exp(-600/T) 

2.2 x 10’ 13 exp(600/T) 

2.33 x I0’"exp(-2814/T) 

1.7 x 10’"exp(-1800/T) 

1.77 x 10’"exp(-2890/T) 

2.14 x I0' l2 exp(-1408/T) 

3.7 x 1 O' 12 exp(240/T) 

1.8 x I0’"exp(110/T) 

6.5 x 10 l2 exp(120/T) 

1.2 x 10 l3 exp(-2450/T) 

1.4 x 10’ 10 

O, 1.91 x I0‘ ,3 exp(-1696/T) 

1.0 x 10" 

2.3 x 10 " 

2.0 x I0’"exp(-3000/T) 

2.0 x ]0'"exp(-3700/T) 

1.8 x 10‘" exp(-390/T) 

3.0 x 10’ 17 

4.02 x I0' l4 exp(317.7/T) 

1.55 x 1 O' 13 exp(-2288/T) 

4.3 x 10' 12 exp(-l 148/T) 

8.59 x 10 " 

1.4 x 10’" 

3.0 x 10' 1J exp(-7000/T) 

3.17 x 10' l2 exp(-4455/T) 

1.2 xlO’ 15 

1.23 x I0‘ l2 exp(-316.8/T) 

4.44 x 10'“ T° ** exp(94/T) 


Technology (NIST) or are extensions to higher or lower tem- 
peratures. Rate constants for the reaction of HNO, with H and 
O (reactions (28) and (29)) are taken from data by Tsang el al 
[1991], 

Reduced Reaction Mechanism 

In the intermediate regime of the wake the jet plume was 
entrained into the wing tip vortex. To model the impact of the 
wing vortex on ozone in this region, large eddy simulations 
(LES) may be required in order to model the vortex properly. 
Solving the flow field by LES would demand significant com- 
puter time, making it impossible to incorporate a detailed 
chemical mechanism. A reduced mechanism capable of de- 
scribing the key chemical reactions is necessary. A reduced 
chemistry is feasible if certain intermediate species reach the 
quasi-steady state when their production rates are nearly equal 
to their consumption rates. Consequently, the concentration of 
a quasi-steady state can be expressed in terms of major spe- 
cies. We identify a species in the quasi-steady state when it 
satisfies the criterion 


* max(|a>f |,|<^ |) 


(14) 


where d>[ and d>l are the production and consumption rates 
of species k, and Xt is the mole fraction of species k. Because 
of the rapid temperature changes in the near field of jet plume, 
a species that does not satisfy the steady state criterion in the 
near field may reach the steady state in the far field. One may 
construct different reduced mechanisms suitable for different 
jet regimes. 

Development of reduced chemistry has been carried out 
using the PSR with a computer-assisted reduced mechanism 
code [He/r el ai., 1993] for the near field of the jet plume 
where large temperature changes take place Assuming quasi- 
steady state for H, HOj, HA, SO, H 2 , HSO,, NO,, N 2 0 5 , and 
O (9 species out of a total of 24 species), we obtained the fol- 
lowing eight-step reduced mechanism, 

(G-Rl) 0, + S0 2 =0j + S0, 

(G-R2) O, + SO, = 20, + SO, 

(G-R3) O, + 20H = H 2 0 + O, 

(G-R4) NO + O, = O, + NO, 

(G-R5) O, + NO + HNO, + SO, = OH + 2NO, + SO, 

(G-R6) 0 2 + S0 2 + HN0 2 = 0H + N0 2 + S0, 

(G-R7) H 2 0 + S0 3 =H 2 S0 4 
(G-R8) O, + CO = O, + CO, 


The rates for the global reactions are expressed in terms of the 
elementary reaction rates as 

<®G-R1 = ®1 + a>S+(»7+<D|| + <B|,+ <B 1 4+0)|j + ©2, + ©j, + 0) 2 s 

— ©33 + ©37 + CD42 — 0043 + ©49 + © 

©G-R2= ©3“ ©4 + ©g+ ©10“ ©13“ ©, 6 - ©20+ ©25 + ©27 

+ ©29 + ©30+ ©32“ ©42 “ ©43 “ ©44 + ©4g ~ ©50 ~ ©51 

- ©52 - ©53 + ©*5 + ©57 - ©j, - ©go 

©G-R3“ ©6+ ©9+ ©10+ ©11 + ©12+ ©17+ ©30+ ©32+ ©47 
©G*R4“ ©19+ ©20 “ ©22 + ©23 “ ©24 “ 2 ©23 - © M - © 27 - ©^ 

+ ©50 + ©51 + ©52 + ©53 “ ©54 “ ©J7 
®OR5 = ©21 “ ©23 + ©25+ ©24+ ©27“ ©52“ ©53+ ©57 
©G-R4“ ©29+ ©29+ ©30“ ©51 + ©54 
©G-R7 = ©39 
©G-RI “ ©41 + ©60 
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Table 2b. Gas Phase Reaction Mechanism: Three-Body Combination and Thermal 
Decombination Reactions 


Reaction 


Low-Pressure Limit ko 


High-Pressure Limit k* 



5.2 x 10* 35 exp(900/T) 

5.51 x 10‘ 31 r 1J exp(l52.6/T) 
1.30 x 10* M T l 
1.05 x 10°°T* 7 expa44.2/T) 
1.28 x 10' a T l ,3 exp(-100.6/T) 
1.54 x 10* M T lll exp(295.3/T) 
6.58 x 10 s T°* 

Oi l .9 x 10* 33 exp(980/T) 

9.09 X 10‘ M r M 


4.2 x I0‘ ,2 exp(-183/T) 


(R50) NO + O — ►NCb 
(R51) NO + OH — ► HNO: 1.91 « 10 J ‘T ,J, e*p(164.8/T) 

(R52) NO, + O — » NOj 6.21 x O' 23 T 5 23 exp(-570.5/T) 

» HNOs 3.50 X 10 “ T 52 exp(-840/T) 


(R53) NO, + OH ■ 
(R54) NO, + NO, 


► N,Oj 


(R55) NtOs ■ » NOi + NO] 

-> OH + NO 
■+ OH + NO, 


(R56) HNO, ■ 
(R57) HNOj 


M 


(R58) SO 2 + O >S(h 

(R59) SO; + OH — HSO, 


6.98 x 10'“ r M 

1.5 x 10 -4 exp(-9104/T) 

5.00 x 10* T 5 * exp(-25340/T) 

1.2 x 10-W-22270/T) 

1.49 x I0'”exp(-601/T) 

1.97 x I0'“exp(867.3/T) 


7.5 x 10'" 

2.9 x I0 '*exp(5821/T) 

2.69 x 10''°exp(-75/T) 

1.5 x 10 " 

2.2 x I0' l> exp(600/T) 

5.43 x 10' 12 T° J 

5.1 x I0' ,2 exp(l59.8/T) 

2.65 x 10'" exp(-25.7&T) 

1 ,84xl0' 12 exp(587.2/T) 
6.40xl0' ,J f : 

5.9 x 10 u exp(-9370/T) 

1.20 x lO 1 *"!” 1 D exp(-25010/T) 
2.14 x 10 l, exp(-24270/T) 

1.19 x 10'* 

2.0 x 10' 12 


(R60) CO + O — — *■ CO, 


2.49 x 10'” exp(-1550iT) 2.66 x 10 l4 exp(-1459/T) 


The effective second-order reaction rate constant is [Demon el al . , 1992] 


k(M,T) = 


k 0 (T)[M] 


\ + (k 0 (T)m/k m (T)) 


:K 


[lo8LO<*o( 7 'XWl/*-(r))] 2 r 


cm 5 s ' 1 and s ', respectively. 


Although the species in the quasi-steady state do not appear 
explicitly in the global reactions, their kinetics effects are in- 
cluded in the elementary reaction rates. The concentrations of 
quasi-steady state species are computed by iteration using the 
following expressions: 

[HS ° 3l_ *40[O 2 ] ’ 


[N°3] = + fc2j [ N o 2 ] + k 27 [OH] + ks4[N0 2 ] + *2i[NO] 

A = Jt 31 (0][HN0 3 l + fc»[N 2 O s ]+ k 52 [0][N0 2 ] 

+ k23[N02][03] + k 32 [OH][HN0 3 ], 


W-f- 


[O] 


k 3 if0 2 lfS01 + * 3 [H1[QH] + fc6[H][HQ 2 ] -t- k|o[OH][OH] 

D = k ti33 {S0 2 ] + Jlr 3 , [SOj] + WHIM + * 42 [0][M] 

+ *i6[H 2 O 2 ] + k 60 [CO] + k 31 [HN0 3 ] + k 2S [HN0 2 ] 

+ k 2 g[N0 3 ]+ti 3 [H0 2 ] + k 5 s[S0 2 ] + ki[0 3 ] + k-,[OH] 

+ k 52 [N0 2 ]+k5o[NO] + knl N0 2 ] + fc 43 [Oj]; 

^6.33[G] [S0 2 1 


[SO] ~ k M [0 3 ] + k 35 [OH] + fc 36 [N0 2 ] + k 33 [0 2 ] ’ 


[tf 2 l = 


B = Jt 35 [OH][SO] + k 7 [OH][0] + &9[OH][H 2 ] 

+ *4i[OH][CO], 

C = ^[H] + k 3 [OH] +k44[OH M l +k i8[H20 2 ] +*29[HN0 2 ] 
+ k 17 [H 2 0 2 ] +k 6 [H0 2 ] + k 5 [H0 2 ] + k 4 7[OH] + k 4 [H0 2 ] 

+ *2(Oj] + ^2 4 (H0 2 ] +fc 4S [0 2 ]; 


MOH] 

E = k^lHKH] + k 3 [H][OH] + ki g [H][H 2 0 2 ] 
+ k 29 [H][HN0 2 ] + k 5 [H][H0 2 ], 


[H 2 0 2 ] = 


fc 49 [H0 2 ][H0 2 1 + fc| 3 [HQ 2 ]fHQ 2 ] + fcutOHHOH] 
k 18 [Hl + k 16 [0] + k n [H] + k 12 [OH] 
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[ho 2 ,4 

F = -t 18 (H][H 2 0 2 ] + ^ l 6 [0][H 2 0 2 ] + it 27 [OH][N03] 

+ ^ 12 [OH][H 2 O 2 ] + ^(OH][O 3 ] + Jt 4 0 [O 2 ][HSO 3 ] 

+ ^45[H][0 2 ], 

G = * 6 (H] + * 5 (H] +MH) + k u [0] + k 49 [H0 2 ] 

+ *1 5 [H0 2 ] + k H [0 3 ] + *n[OH] + Atjo [NO]; 

[N 2 o 5l = MNQ 2 ][No 3 ] 

*55 

Results and Discussions 


action (30) becomes less important for mixtures with high 
free-stream/exhaust ratios. Reaction (19) becomes the domi- 
nant route for NO 2 production for mixtures with large 
free-stream/exhaust ratios as more O 3 is present in the system. 
Reaction (51) has positive sensitivity with respective to NOj 
as it produces HN0 2 which then forms N0 2 via reaction (30). 
Figure 6 shows the sensitivity information of O 3 over all mix- 
tures. As expected, reaction (19) NO + Oj-* NOj + O 2 is the 
most important step in controlling O3. For mixtures with low 
free-stream/exhaust ratios, 0 3 is produced via reaction ( 43 ) O 
+ Oj + M-*Qj+M, and N0 2 is converted to NO via reaction 
(22) NO 2 + O — * NO 4 - O 7 . When free-stream air is increased, 
the importance of these two reactions diminishes due to the 
consumption of the O radical. 


Results are presented for a Mach 2.4 HSCT flight flying at 
18 km and N85. The initial exhaust and ambient conditions 
are taken from Brown et al. [19%]. A jet exit radius of 0.93 m 
[ Kdrcher , 1 994] is used in the simulations. 


Sensitivity Information and Key Reaction Paths 

Sensitivity information is useful for identifying the control- 
ling chemical steps responsible for the creation/destruction of 
a certain species. The sensitivity coefficients are displayed as 


B -A.* 3 ** 


(15) 


where fi ik is the sensitivity for change in the mole fraction of 
the foh species, due to a small change in temperature- 
independent prefactor, A h of the ith reaction rate k, = A ( F 
exp (-EJT). We used the ratio (mole base) between the free- 
stream (ambient air) and the exhaust gas from the engine to 
describe the different mixtures. A large free-stream/exhaust 
ratio corresponds to a mixture highly diluted with ambient air. 
The PSR model has been used to assess the degree of changes 
in ozone and NO* in the engine exhaust gas when it is per- 
fectly mixed with a certain amount of ambient air over resi- 
dence times up to 100 s. Figure 4 presents a set of sample 
sensitivity data of NO versus residence time for four free- 
stream/exhaust ratios at 0.1, 1, 10, and 100. As seen from the 
plots, both reaction (19) NO + 0 3 -► NO: + O? and reaction 
(51) NO + OH + M HN 02 + M are responsible for deplet- 
ing NO. When there is a large amount of air mixed with the 
exhaust gas, reaction (19) becomes the dominant step in con- 
verting NO into NOj. Since reactions (53) OH + N0 2 + M 
HNO 3 + M and (30) HNO 3 + OH —> H 2 0 + NO 3 are compet- 
ing with reaction (51) for OH, increasing their rates will re- 
duce the effectiveness of reaction (51) in converting NO into 
HNO}. Consequently, reactions (53) and (30) have positive 
sensitivity values. However, the importance of reactions (53) 
and (30) diminishes as more air is mixed with the exhaust gas 
since OH concentration becomes small compared to O 3 in 
mixtures with high free-stream/exhaust ratios. 

Next we examine the sensitivity information for NO, as 
presented in Figure 5. As revealed in the sensitivity data, re- 
action (53) OH + NO: + M — * HNO 3 + M is the major chan- 
nel for converting NO: into a more stable species HNO 3 . 
Consequently reaction (53) has a negative sensitivity for 
all mixtures. The sensitivity data also show that NC^ is 
produced via reaction of NO with O 3 through reaction (19) 
and via reaction of OH and HNO : through reaction (30). Re- 


Evaluation of Reduced Chemistry 

First, PSR is used for evaluating performances of the eight- 
step reduced by comparing concentrations of NO, N0 2 , and 
O 3 over a wide range of free-stream/exhaust ratios for resi- 
dence times from 0.01 s to 100 s. In all cases, excellent 
agreement is found between the results obtained with the de- 
tailed mechanism and those with the reduced chemistry. Fig- 
ure 7 presents sample comparisons showing almost identical 
results when reduced chemistry is used; however, the com- 
puting time is decreased by a factor of 3 . 

Second, to assess the transient responses of reduced chem- 
istry to time-varying mixing with ambient air, WMR simula- 
tions were performed with air injection rates estimated on the 
basis of the PDF jet simulations. As described above, the 
WMR simulations are similar to the commonly used La gran- 
gian box model except that the temperature is not prescribed 
but is solved for directly using the energy equation. Figure 8 
presents the time evolution of major nitrogen, sulfur, and oxy- 
gen species for the engine jet exhaust over the first 2 seconds 
after the exhaust nozzle, showing excellent agreement. To 
further evaluate the performance of reduced chemistry with the 
micromixing effect, the PaSR model is computed with identi- 
cal conditions as those in the WMR model, and the results are 
presented in Figure 9. Consistent with the comparisons ob- 
served in the WMR calculations, the performance of the re- 
duced chemistry remains quite satisfactory in spite of the fact 
that a wider range of mixtures have been computed since the 
reactor is not well mixed 

Third, Monte Carlo PDF jet simulations have been per- 
formed and Figure 10 presents the predicted centerline tem- 
perature evolution obtained with the effects of microscale 
mixing. The computed temperature is seen to decreases rap- 
idly in the near field, manifesting fast entrainment of cool am- 
bient air. Figure 1 1 presents the centerline mole fractions of 
major nitrogen, sulfur, and oxygen species evolution obtained 
with the detailed mechanism and the eight-step reduced 
chemistry. It shows good agreement except for O 3 at the early 
jet stage. This is due to the steady state assumption for the O 
radical which does not achieve steady state at early jet stage 
where it contributes to the production of (X through O + O? + 
M — > O 3 + M. 

Impact of Microscale Mixing on Chemical Kinetics 

A parametric study of the influence of various mixing proc- 
esses on predicted chemistry has been carried out with the 
models listed in Table 1 . To quantify the impact of different 
mixing processes on NO* , HNO 3 , and O 3 in the jet plume, we 
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Figure 4. PSR model predictions of sensitivity of NO species to elementary steps at free-stream and exhaust 
mole ratios of 0.1, 1, 10, and 100. 


examine the species response normalized by NO, flux at the 
nozzle exit. Because of the small amount of NO, entrained 
into the exhaust plume the total NO, mass flux remains nearly 
constant. The normalized species response is defined as 

normalized species response > \ l °) 

m nox 

where m(t) is the mass flux of species / evaluated over the 
cross section of the jet at a certain downstream location, and 


HH total is ^ sum to® flux of a ccrtaui splits origi- 
nated from the nozzle exit and the mass flux entrained into the 
jet, that is, A species re- 

sponse value of zero indicates that the total amount of species 
is conserved Negative or positive species responses corre- 
spond, respectively, to destruction or creation of species due 
to chemical reactions. According to convention, the mass of 
NO, is reported as the equivalent mass of NO?, which is ex- 
pressed as rrtytck * m NOj + ^ NO] / W No m No » where W is the mo- 
lecular weight of a species. 
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Figure 5 . PSR model predictions of sensitivity of NO2 species to elementary steps at free-stream and exhaust 
ratios of 0. 1, 1, 10, and 100. 


ics ongi- 
’i into the 
ecies re- 
f species 
2s corre- 
xies due 
mass of 
ch is ex- 
the mo- 


Figure 12 presents the predicted NO,, HN 0 3 , and 0 3 re- 
sponses obtained with various numerical models. As seen 
from the NO, and HN 0 3 responses, rapid conversion of NO 
and N 0 2 into HNO3 occurs in the near field due to reactions of 
OH + NO z + M -► HN 0 3 + M and NO + 0 3 -► N 0 2 + 0 2 . 
After 0.5 s, the NO, response remains nearly constant, indi- 
cating negligible changes in total NO, due to chemical reac- 
tions. The HN03 response is directly correlated with the NO, 
response as its production is via the reaction between OH and 
NO2. The 0 3 response shows that a small amount of 0 3 is 


produced via reaction O + Oj + M — ► 03 + M at early jet stage, 
leading to a peak O3 response at about 0.05 s (-70 m from 
nozzle exit plane). After the peak, entrainment process be- 
comes important, and the 03 response starts to decrease due to 
the dominant destruction reaction NO + 03 -> NO 3 + O: as 
more 0 3 is entrained. The predicted 0 3 depletion per unit mass 
of NO, increases roughly linearly with time or roughly with 
downstream distance. The constant destruction of 0 3 per unit 
mass of NO, exhaust is expected as NO, plays the catalytic 
role in consuming 03. 




12,880 


WANG 


AND CHEN: MODELING OF TURBULENCE AND CHEMISTRY IN PLUME 





Overall, conversion of NO x into HNO 3 decreases when 
more detailed features of mixing processes are included in the 
model. The WMR model (similar to a Lagrangian box model) 
gives the largest NO, conversion into HNO 3 since both mac- 
roscale and microscale mixing processes are assumed infi- 
nitely fast When microscale mixing is accounted for in the 
PaSR, the NO, response reduces by 3.5%. The PDF jet 
simulation with well-mixed assumption for each grid cell 
(similar to the SPF-II code) includes the macroscale mixing 
process that brings ambient fluids into the jet plume through 


turbulent convection and stirring. The differences in the NO, 
responses between the WMR and the PDF jet model with 
completely mixed cells are about 25%. When microscale 
mixing is included by the modified Curl’s mixing model, the 
NO, response reduces further by 5%. Consequently, the 
amount of HNO, produced is largest when infinitely fast 
mixing is assumed. The influence of mixing processes at dif- 
ferent scales on the predicted O 3 depletion is seen similar to 
the NO x response except that the trend is reversed. Within the 
first 0.05 s, O 3 is produced by O + 02 + M -► 0 3 + M. The 
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Figure 7. Companson of PSR model predictions of NO, N0 2l 
and O} mole fractions with detailed and reduced mechanisms, 
showing good agreement over a wide range of mixtures of ex- 
haust and ambient air as a function of residence time. 



tom* (mc) 


Figure 9. Comparison of partially stirred reactor (PaSR) pre- 
dictions of major nitrogen, sulfur, and oxygen species evolu- 
tion with detailed and reduced mechanisms, showing good 
agreement 


timescale of chemistry for the radical O is comparable with 
that of microscale mixing. A strong coupling between reaction 
and mixing is noticed as the amount of O* produced increases 
with mixedness of the jet plume. After about 0.05 s, the de- 
struction reaction N0 2 + Oj — > NO + 0 2 becomes dominant 
due to fast mixing of jet plume with ambient air. Since the 
amount of O 3 depleted via this reaction is almost the same, the 
O 3 response curves for each model evolves in an almost par- 
allel fashion. 

A comparison of the predicted centerline mole fractions of 
major species and the relative differences obtained using the 
PDF model with/without Curl’s mixing model is presented in 
Figure 13. The comparison shows the influence on chemical 
kinetics due to microscale mixing only. When the fluids are 
assumed well mixed at the microscale, the point at which en- 
trainment becomes important starts earlier. The impact of mi- 
croscale mixing at axis is noticeable after 50 m downstream, 
where the timescale of chemical reactions for radicals, such as 
OH and HOj, is comparable or faster compared with that of 
microscale mixing. The impact of micromixing on radicals 
affects concentrations of other species since hydroxyl radical 
is the driver for major chemical reactions in the near field. In- 
clusion of microscale mixing increases predicted H 2 S0 4 before 
the axial distance of 500 m and the increase reaches about 
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Figure 8. Comparison of WMR predictions of major nitrogen, 
sulfur, and oxygen species evolution with detailed and re- 
duced mechanisms, showing good agreement. 
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Figure 10. Centerline temperature evolution for probability 
density function (PDF) jet simulation. 
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Figure 13. (a) Centerline mole fractions of major nitrogen, 
sulfur, and oxygen species for PDF jet simulation with differ- 
ent treatment for microscale mixing, (b) Relative difference of 
major species between the two models Relative difference 
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Figure 12. NO., HNOj, and Oj responses for different nu- 
merical models to explore the macromixing and micromixing 
effects. 


20% at 100 m. This may have a potential impact on predic- 
tions of sulfate aerosol formation which is a subject to be 
studied in the future. 

Conclusions 

Using a PSR model with a detailed reaction mechanism, 
the key reaction steps controlling NO,, O 3 , and OH evolution 
have been identified in engine exhaust gas mixed with various 
amounts of ambient air in the stratosphere. On the basis of the 
detailed mechanism, an eight-step reduced chemistry has been 
developed and tested extensively using PSR, WMR, PaSR, 
and the PDF jet simulation. The performances of the reduced 
chemistry are found quite good over a wide range of mixing 
conditions. 

The relative importance of turbulent mixing at macroscale 
and microscale on chemistry is investigated by a parametric 
study using various models for the near-field jet plume. For a 
typical HSCT flight condition, inclusion of macroscale mixing 
leads to about 25% lower NO, conversion to HNO 3 than that 
predicted by the WMR model (similar to a Lagrangian 
model). In regions close to nozzle exit, O 3 is produced via the 
reaction O + O, + M -► O 3 + M. The amount of 0 3 produced 
increases with the degree of mixedness inside the exhaust 
plume. As such, the WMR model yields highest 0 3 peak re- 
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sponse in the near field in comparison with other models. The 
impact of macroscale mixing on chemical kinetics is found 
most significant in the early development of exhaust plume 
within 0.2 s. Inclusion of microscale mixing causes a slight 
delay in the mixing of exhaust species with ambient air. This 
impact results in about 5% changes in the overall species re- 
sponses of O}, NO„ and HN0 3 across the jet However, the 
influence of microscale mixing on local species concentration 
is found to cause nearly 50% differences for highly reactive 
radicals in the near field. With the microscale mixing included 
in the numerical model, the predicted centerline peak H 2 S0 4 
concentration increases at locations where sulfate aerosol for- 
mation could be highly active. The potential impact of micro- 
scale mixing on aerosol formation is a subject of future study. 
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Numerical Studies of Near-Field 
Plume- Vortex Interactions 


Junxiao Wu* and Suresh Menon* 

School of Aerospace Engineering 
Georgia Institute of Technology 
Atlanta , Georgia 30332-0150 

In this study, fluid mixing and chemical reactions in the near-field wake of an aircraft 
has been investigated. The plume-vortex regime is studied using a 3D parallel LES code. 
Simulations of engine plume interaction with the wing-tip vortex of B737 and B747 have 
been carried out and some comparisons with data have been used to demonstrate the abil- 
ity of the LES method. Inclusion of detailed gas phase kinetics and simple heterogeneous 
processes have been carried out to evaluate the effect of turbulent mixing on the chemical 
processes in the plume-vortex regime. Finally, to determine the impact of temporal as- 
sumptions (often used in such flow simulations) comparison with spatial simulations are 
carried out. Many features of the computed flow field show reasonable agreement with 
experimental data. The entrainment of the engine plumes into the wing-tip vortex and 
also the partial detrainment of the plume were numerically captured and analyzed. The 
impact of fluid mechanics on the chemical process has been estimated. Results indicate 
that a significant difference between spatial and temporal simulation exists (especially in 
the predicted SO z concentration). This has important implications for the prediction of 
sulfuric acid aerosols in the wake and may partly explain the discrepancy between past 
numerical studies and the data. Analysis of the results also show that spatial simulation is 
more suitable for the near-field interaction process. For example, it was determined that 
detrainment (observed in the B747 contrail) was only captured using the spatial model. 


1 Introduction 

The wing-tip vortices shed from an aircraft's wings 
are known to pose a danger for aircrafts in the wake 
and are considered a major hazard near airports. An- 
other consequence of the shed vortices that is perhaps 
less appreciated is their ability to entrain the hot ex- 
haust of the engine during cruise flight in the upper 
atmosphere. The resulting species-rich and relatively 
long-lived contrails have been of increasing concern 
and a focus of research in recent years because the 
species entrained into the vortex wake can react with 
the ambient species thereby changing the local chemi- 
cal balance. These studies have been motivated by the 
projected increase in future subsonic flights in the up- 
per troposphere by the proposed Advanced Subsonic 
Transports (AST) operating in the transatlantic and 
Pacific rim flight corridors and by the possibility of fu- 
ture supersonic flights in the lower stratosphere by the 
High Speed Civilian Transport (HSCT). 

The near-field region of the aircraft's wake can be 
divided roughly into two parts: a plume regime and a 
plume-vortex regime. The plume regime is identified 
as the initial region behind the aircraft where the en- 
gine exhaust plume grows and mixes with the ambient 
air. This regime is typically very short, approximately 
l sec after exhaust or around 200 m downstream of the 
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engine exit plane. Further downstream is the plume- 
vortex regime where the exhaust plume is entrained 
into the rolled-up wing tip vortices. This regime is 
around 1-2 km long and approximately 10 sec down- 
stream of the exhaust. Due to the relatively high tem- 
perature of the engine exhaust and the species content, 
significant reactions can take place in the plume- vortex 
regime. The region downstream of the plume- vortex 
region is the wake region where the vortex wake breaks 
up due to (Crow) instability and mixes with the ambi- 
ent air. This region is very large extending over 20 km 
and the mixing process can t ake many days depending 
upon the ambient conditions. 

In the present study, the near-field region, especially 
the plume-vortex regime is of interest. In particular, 
we are interested in the chemical processes that can 
occur in this region due to the reaction of the emit- 
ted species with the ambient air (in particular, ozone). 
Key species that are of interest here are nitrogen ox- 
ides (NO x ) (generated from the oxidation of molecular 
nitrogen at high temperature), carbon monoxide(CO) 
(due to incomplete combustion), methane (C H 4), non- 
methane hydrocarbons(//C) and soot particles. The 
role of these trace species in the atmospheric pho- 
tochemical and radiative processes is not very well 
understood. For example, the concentration of tro- 
pospheric ozone, an important greenhouse gas, is de- 
pendent on transport across tropopause and by chem- 
ical production and destruction involving reactions of 
NO x i non-methane hydrocarbons, carbon monoxide. 
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and hydrogen oxides ( HO r )• Furthermore, the radia- 
tive balance of the atmosphere can be affected by the 
accumulation of carbon dioxide and water ( CO 2 and 
H 2 0 are the primary species in the engine exhaust) 
and by the aircraft-induced changes in aerosols and 
clouds. In addition, the combustion aerosols (soot) 
and the binary HnSOij HoO droplets in the plume can 
be efficient cloud condensation nuclei and serve as sites 
for heterogeneous reactions. 

The ozone balance in the upper troposphere and 
lower stratosphere can be perturbed by an increase in 
XO r and by heterogeneous reactions that take place 
on the surface of aerosols . 1 ’ 2 In the troposphere, N O x 
acts as a catalyst to produce O3 in reactions involv- 
ing the oxidation of methane and other hydrocarbons. 
However, in the stratosphere, where O3 is a barrier 
to ultraviolet light, \O r catalytically destroys O3 in 
react ions involving 0 atoms. Thus, the impact of emit- 
ted species on ozone is a function of many variables. 
This can complicate long-term prediction of upper at- 
mosphere chemical balance. 

A variety of models have been developed to assess 
the effect of aircraft flights on the atmosphere. These 
models generally fall into two categories: " small-scale" 
and "global" models. "Small-scale" models treat the 
near-field regime of the exhaust wake by ignoring pro- 
cesses that operate on larger time and space scales. 
Global models . 3 on the other hand, attempt to repre- 
sent the large spatial scales and ignore, or represent 
verv crudely, the small scales processes. To obtain 
some sort of continuity, some of the global models 3 use 
the engine exhaust conditions or the data at the end of 
the plume-vortex regime as input to their large-scale 
models. However, as shown recently , 4-6 the concentra- 
tion of species (e.g., chemical radicals OH , HOn and 
soot) can undergo very rapid changes in the near-field. 
This can significantly modify the input conditions for 
the far field simulations. Thus, to obtain accurate pre- 
diction of the global chemical balance due to aircraft 
emission, the chemical processes in the plume- vortex 
regime is very important. 

Models for the near-field regime have been de- 
veloped by many researchers. The simplest model, 
often called the "trajectory box" model , 78 ignores 
flow transport and focuses primarily on the chemical 
species evolution. This type of model has been very 
popular since detailed chemical processes can be stud- 
ied due to its computational efficiency. Models that in- 
corporate both chemical and transport processes have 
also been developed. However, some approximations 
were employed to reduce the computational expense. 
For example, some models 4,9,10 used the well known 
SPF-II code 11 which is an efficient parabolic Navier- 
Stokes model. Gas phase kinetics, binary H 2 SO 4 - 
H 2 0 nucleation and aerosol coagulation have been 
studied using this model. Gamier et al . 12 used an 
integral model and a two-dimensional direct simula- 


tion of convection-diffusion equation to investigate the 
evolution of the mixing along the wake of a typical 
large transport aircraft. Menon and \\u 6 developed 
a new mixing model that included both micro- and 
macro-scale mixing effects and Wang and Chen 0 used 
a scalar probability density function model which in- 
cluded the large-scale mixing effect but ignored the 
micro-scale mixing effects. Comparison of these mod- 
els was recently carried out 6 and it was demonstrated 
that inclusion of both macro- and micro-scale mixing 
results in 30% less 0 3 and 15% less \ 0 2 depletion in 
the near field regime . 6 

The hydrodynamics of the plume-vortex regime is 
also quite complex. For example, the competition 
between the vortex centripetal force and gravitation 
buoyancy of hot plume can lead to a different final al- 
titude of deposition of emissions . 9 This implies that 
evolution of the plume occurs in a highly unsteady 
manner. However, many of the models noted above 4 0 
are steady-state models and therefore, are unable to 
account for the unsteady mixing effects. 

To study unsteady jet-vortex interactions, unsteady 
methods have also been used. For example, direct nu- 
merical simulation 13 and large eddy simulation (LES) 
of trailing vortices 14 showed that t here exists strong 3- 
D instability in vortex wakes with an axial deficit. The 
vortex instability results in vortex-stretching which 
may have a substantial effect on the mixing of jet 
plume with the ambient air. Gerz and Ehret u used 
LES to simulate the vortex roll-up process of the wake 
of a B-747 one second after exit. They showed that 
temperature, relative humidity peak at the center of 
the wingtip vortex, the wing boundary layer turbu- 
lence and the ambient turbulence all can affect t he 
plume-vortex interaction process. 

Unsteady simulations are computationally very ex- 
pensive and therefore, most past studies employed 
some simplifying assumptions. For example, only lim- 
ited chemical effects have been investigated so far lh 
and all simulations reported so far employed periodic 
boundary conditions in the streamwise directions (i.e. 
simulated temporal evolution instead of the actual 
spatial evolution). Temporal approximation is truly 
justified only far downstream where the interaction 
between the jet plume and wingtip vortex has been 
completed. It has not yet been shown what features 
of the flow and the accompanying chemical processes 
are different and/or modified by this temporal assump- 
tion. This issue is important since it has been shown 
that the entrainment predicted by temporal and spa- 
tial simulations differ markedly due to the asymmetry 
between the high speed and the low speed sides of the 
shear layer . 17 The impact of this asymmetric entrain- 
ment on the chemical processes in the plume- vortex 
interaction regime has not been addressed so far. 

This paper will address many of these issues and will 
attempt to quantify the impact of simulation method- 
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ology (i.e.. spatial or temporal) on the prediction of 
the species distribution. Furthermore, the effect of 
the fluid dynamic interaction process on the chemical 
kinetics and on the concentration of species such as 
.VO, S0 3 and 0 3 in the plume-vortex regime will be 
addressed in the presence of both gas and heteroge- 
neous kinetics. 


2 Simulation Model 

The simulation model used in the present study em- 
ploys the LES method. In LES, scales smaller than 
the grid size 6. v are eliminated by applying an ap- 
propriately chosen low-pass filter G to the flow field 
variables. Any filtered field of a quantity / is given 
as: 7 = / d G(x - x f )f(x')dx'. A box filter is used 
since it is appropriate for finite- volume schemes. 13 Us- 
ing the above definition, the original field / can be 
decomposed as / = 7 + /'• For the compressible 
flow equations, a Favre-filter is more appropriate. 19 
In this case, the Favre-filtered variable is defined as: 
/ = £l. The details of the filtering operation on the 
Navier-Stokes equations are given in numerous pub- 
lications 20,21 and therefore, avoided here for brevity. 
The LES equations for conservation of mass, momen- 
tum. energy and species are, respectively: 
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In the above equations, p is the density, u, is the 
velocity component in ith direction, p is pressure, T is 
temperature and k is the thermal conductivity. Also, 
Y k . ^j k< Vf are, respectively, the kth species mass frac- 
tion, mass reaction rate and diffusion velocity. The 
kth species diffusion velocity is determined using Fick's 


Law of diffusion as: \j k = -Djf-JYk. " here D is the 
mixture diffusion coefficient. Also. r~ is the viscous 
shear stress determined in terms of the resolved veloc- 
ity gradients. 

The filtered total energy E per unit mass is defined 

83 „ 1 

E = e 4- -itiiii + k sg$ (5) 

where the subgrid kinetic energy, k sgs , is defined as 

k’ 33 = ^(upUi - «>'.) ( ( >) 

Finally, the LES equation of the state equation is 
p = pRT (7) 


In the above equation, R is the mixture gas con- 
stant given as YkRu/Wk where R u and \\ k are 

the universal gas constant and the species molecular 
weight for kth species respectively, and A is the total 
number of species. 

2.1 Subgrid Closure of the LES Equations 

The above LES equations contain many terms that 
arise due to the filtering process and represent the 
effect of the unresolved subgrid fluctuation correla- 
tions on the resolved motion. Models for these un- 
known terms must be devised to complete the LES 
formulation. The key subgrid term that requires 
determination is the subgrid Reynolds stress term: 
r 3gs = p[EiUj - Uiiij] and the subgrid enthalpy flux 

Hi sgs = p[Hui - Hui]. Here. H = E + p/p is the 
stagnation enthalpy per unit mass. 

A subgrid eddy viscosity closure is typically chosen 
for closure of the momentum and energy subgrid fluxes 
since the small scales primarily provide dissipation for 
the energy transferred from the large scales. To de- 
termine the eddy viscosity, characteristic length and 
velocity scales must be defined. The grid size is typi- 
cally chosen as the length scale and the velocity scale 
in the present study is obtained by solving a transport 
equation for the subgrid turbulent kinetic energy. k* ,J9 . 
Without going into the details, 21,22 the final equation 
for subgrid kinetic energy is given as: 


3pk sgs 
dt + 


dx x 


d dk $gs 

— — ) 
dx t Ox; 

+ Pk9fJ9 ~~ Dktgf {$) 


In this equation, the three terms on the right-hand-side 
are respectively, the transport term, the production 
term, Pksgs = -r/f* diti/dxj and the dissipation term 
Dksgs — c £ p(k S93 ) 3 / 2 / A. Here, A is a characteristic 
grid size. Once the subgrid kinetic energy is deter- 
mined the subgrid stress tensor, r* gs can be obtained 
as: 

r*f* « -2piy t (S,j - jSfcfciij) + ^pk ,y ’Sjj (9) 
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where 5 0 = 5 (ft + §£) the resolved rate-of-strain 
tensor and i/ f is the subg rid e ddy viscosity which is 
approximated as v t % c u \Jk s $ s A, where c v is a model 
coefficient. Similarly, the subgrid enthalpy flux term is 
modeled as: pHiii — pHiii = — c e py/k*9* A|^. Here, 
c u , c ( and c € are unknown coefficients that must be 
specified to close the system of equations. Typically, 
c e is modeled as c^/Pvt where the turbulent Prandtl 
number was chosen as 0.9. The method used to eval- 
uate these coefficients is described in the next section. 

Closure of the s pecies equation requires closure of 
the subgrid term: p^kUj — P^k l *j * A subgrid diflfusiv- 
ity model (as the eddy viscosity model given above) 
is not considered an excellent choice since for com- 
bustion to occur species must first mix at the small 
scales and then, molecularly diffuse. These processes 
occur at the unresolved scales and therefore, a global 
eddy diffusivity model for the subgrid species flux that 
ignores the subgrid processes requires some justifica- 
tions and can be used only under certain conditions. 
In the present study, the mixing process is dominated 
by the large-scale transport of the plume species into 
the wing-tip vortices. Furthermore, the time-scales 
for the subsequent chemical reactions are quite large 
(i.e., slow chemistry). Thus, the fast time small-scale 
mixing effects can be considered completed before the 
chemical reactions occur. For these reasons, an eddy 
diffusivity model for the species subgrid flux is con- 
sidered a reasonable first-order approximation. With 
these comments in mind, the subgrid species flux is 
approximated as: 

— — ,, m 
= < 10 > 


where Sc T is the turbulent Schmidt number set to 
unity for the present study. 

2.2 Dynamic Evaluation of the Coefficients 

The model coefficients c*,, c f and c e are obtained us- 
ing a dynamic approach. The details of the dynamic 
approach has been reported elsewhere 21 and therefore, 
is only briefly summarized here. To obtain the model 
coefficients, we assume that the subgrid stresses at the 
grid filter level and the resolved Leonard's stress at 
the test filter level with a characteristic grid size A 
(at twice the grid filter) are similar. Scale similarity 
between r t j and Lij has been observed in some ex- 
perimental data. 23 Using this similarity, the model 
coefficients are obtained from the following relations: 
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The ability of this type of closure to obtain accurately 
the LES field for a wide range of test problems in 
both incompressible and compressible flows has been 
reported in the past. 20,21 It has also been shown in 
these earlier studies that predicted subgrid stresses us- 
ing this model satisfies all the realizability conditions 
proposed by Schumann. 24 This implies that since 
is a positive semi-definite tensor, the modeled term 
for r-f 4 5 should also be positive semi-definite. This re- 
quirement results in the following inequalities: 


r*r > 0, / E {1,2.3}. 

| r ^|2 < r ,gs r jg* JJe {1.2,3}. (18) 

det(r?f *) > 0. 

that must be satisfied. The above inequalities yield 
certain bounds for c u which must be imposed in the 
simulations. The present implementation of the dy- 
namic model ensures that these requirements are satis- 
fied during the entire simulation. Analysis shows that 
this criteria is automatically satisfied in over 95 ( /c of 
the grid points during the entire simulation. 


3 Numerical Method 

An explicit finite- volume scheme that is fourth order 
accurate in space and second-order accurate in time 
was used in the present study. The details of the nu- 
merical scheme has been reported elsewhere 21 *" and 
therefore, avoided here for brevity. The current ap- 
proach is similar to the implementation described by 
Garnet and Estivalezes 25 except that the fluxes are 
modified so that a truly fourth-order accuracy in space 
is obtained at least on an uniform grid. The details of 
this implementation are given in. 22 
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Table 1 Aircraft parameters for the near field 
plume-vortex calculations of B74T and B737 data. 


type 

wing 

span 

maximum 

circulation 

vortex 
core radius 

vortex 

separation 

engines 

location 

B747 

60 m 

600m* /s 

4 m 

54 m 

Dm/ 18m 

B737 

30m 

100 m“/s 

2 m 

26m 

5.35 m 


Table 2 Initial conditions for the near field plume- 
vortex calculations of B747. 



B 

747 


Plume Condition 

Ambient Condition 

Temperature 

500K 

217K 

velocity (m/sec) 

247 

0.0 

pressure (atm) 
mole fraction 

0.17 

0.17 

CO 

2.4xl0 -5 

2.0xl0 -s 

co 2 

3.2xl0~ 2 

3.5x10—* 

Ho 

l.OxlO -6 

5.2x10“ ‘ 

Oo 

0.159 

0.2098 

OH 

l.OxlO -5 

2.8xl0 -13 

HoO 

3.0xl0 -2 

4.8xl0 -6 

No 

0.779 

0.789 

NO 

4.3xl0 -5 

5.9x10"“ 

,\ Oo 

4.8xl0 -6 

4.3xlO~ 10 

N 0 O 5 

0.0 

0.0 

HNO 3 

0.0 

3.3xl0 -s * 

HONO 

0.0 

0.0 

SOo 

6.9xl0 -6 

0.0 

so 3 

0.0 

0.0 

0 3 

0.0 

2.5x10" 7 


The present code has been parallelized using the 
standard Message Passing Interface (MPI) library. 
This allows simulation on a variety of platforms with- 
out any changes to the code. Domain decomposition 
method is employed whereby the computational do- 
main is decomposed into subdomains and assigned to 
different processors which have an identical copy of 
the numerical algorithm. For the current scheme, a 
five point stencil is used and therefore, two rows of 
variables must be exchanged in either direction. Tim- 
ings analysis of this code on various parallel machines 
(e.g., Cray T3E, SGI Origin 2000) shows that the code 
achieves nearly linear speedup as the number of pro- 
cessors is doubled. However, depending upon the total 
number of grid points employed there is an optimum 
number of processors beyond which the message pass- 
ing overhead overtakes the computing time in each 
processor s domain. 

For a 121 x 97 x 97 computational grid a typical sim- 
ulation (using 15 reactive species) on the Cray T3E 
takes around 5000 single processor hours to obtain 
data sufficient for statistical analysis (5 flow through 
time). The CPU time per time step per grid point was 
around 3E-6 sec on 128 processor T3E (or 3.6 sec per 
time step). 

Two types of problems are studied in the present 
study. Both are in the plume- vortex regime. Figures 
la and lb show respectively, the plume- vortex domain 
behind a Boeing 737 and a Boeing 747 aircraft. The 
primary differences are in the span, the number of en- 
gine plumes and the exhaust conditions. These two 
configurations w r ere chosen primarily because there is 
some data available for comparison. 

For slip or symmetry conditions, values are pre- 
scribed for ghost cells such that gradients normal to 
the face at the bounding face are zero. Periodic bound- 
ary conditions are prescribed in the streamwise direc- 
tion (for the temporal simulations). 

For the spatial cases, the inflow and outflow bound- 
aries are more complicated and must be carefully im- 
plemented. We employ a method that attempts to 
provide reasonably accurate inflow turbulence statis- 
tics. 2226 The inflow is critical for the spatial case since 
it determines the turbulence entering the computa- 
tional domain. A turbulent inflow field is generated 
based on a prescribed energy spectrum and turbu- 
lence intensity and then superimposed on the mean 
field. The turbulence spectrum chosen is of the form 
E(jfc) = Ck A exp(-2(k/k p ) 2 ) used by Lee et al 27 where 
constant C* is chosen as turbulence intensity level. The 
turbulence is divergence free and non-periodic in time. 
The field is stored as 3D array from which a 2-D plane 
is introduced at the inflow at every time step. At least 
10 different arrays of 97 x 97 x 97 are stored and to 
ensure that the inflow field is not correlated, at ran- 
domly chosen instants, the phase in the Fourier modes 
were randomized. Additional details of the turbulence 


inflow are given elsewhere.’” 

The outflow' boundary is non-reflecting character- 
istic boundary conditions. 26 The amplitudes of the 
outgoing waves are computed based on the interior 
points and the amplitude of the incoming wave is com- 
puted by specifying the background pressure. Details 
of this type of outflow condition is given elsewhere. 

The wing-tip vortex is modeled by the Lamb-Oseen 
vortex 28 with an axial velocity deficit. This model 
has been extensively used 12 " 1,4 29-31 for earlier studies. 
The location, core radius, maximum circulation data 
of the wing-tip vortex and the engine location and exit 
conditions (Table 1 and 2 summarize these conditions) 
are chosen according to the aircraft configurations. 

4 Results and Discussion 

In this section, we summarize the results obtained in 
the present study. There are three primary objectives 
of this study. First, comparison with available data 
is carried out to demonstrate t he ability of the sim- 
ulation model to accurately capture the evolution of 
the plume- vortex regime. Second, the results obtained 
using temporal and spatial methods are compared to 
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identify what features are not predicted accurately us- 
ing temporal simulations (as noted earlier, nearly all 
previously reported 3D studies employed periodic con- 
ditions in the streamwise direction). The final objec- 
tive is to determine the impact of turbulence transport 
on the formation of SO x and the destruction of ozone. 

Before discussing the details of the plume-vortex in- 
teraction, the effect of different inflow turbulence on 
the flow field is briefly summarized. Figure 2 shows 
the effect of inflow turbulence on the scalar (here SO 2 
concentration is used as a scalar marker) decay with 
increasing plume age in a B747 plume- vortex region. 
Three different inflow turbulence conditions are used. 
Two conditions employed the same energy spectrum 
but with 5% and 10% turbulent intensity, respectively 
while the third case used the same initial isotropic 
spectrum as other two but is introduced into the flow 
field after evolving for certain time with a final tur- 
bulent intensity 0 %. Thus, in the third case, the in- 
flow has evolved to a more realistic turbulent isotropic 
state, whereas in the first two cases the field is isotropic 
but is not realistic turbulence. 

It can be seen that the cases with the same tur- 
bulent intensity yield almost identical results. This 
means that the effect of evolving of the inflow turbu- 
lent field is small at least for the global measure used 
here. (Note that such global measures are of interest 
here). However, it can be seen that turbulent inten- 
sity does have an effect on the flow field. This implies 
that the inflow turbulent intensity must be carefully 
chosen. The turbulence intensity in the vortex and jet 
plume varies with flight conditions such as climate, al- 
titude and ambient flow field. A turbulent intensity 
between 0 % and 10% is very common in the typical 
flight path of interest. Thus, simulations with 5% and 
10% inflow turbulence may be used to quantify an up- 
per and lower bound for the species distributions in 
actual flights. Interestingly, the descending vortex ve- 
locity computed from these three cases showed similar 
results. However, this is understandable because the 
descending vortex velocity is mainly determined by the 
overall vortex circulation and small scale turbulence 
exerts negligible effect on the circulation. The numer- 
ical simulation of q-vortex by Ragab and Sreedhar 29 
also demonstrated this feature. 

4.1 Plume- Vortex Interaction behind a 
Boeing 737 

For quite some time under the Atmospheric Effects 
of Aviation Program, NASA has been conducting mea- 
surements in the wake of actual aircrafts. Some data 
based on ground based LIDAR measurements in the 
wake of a Boeing 737 is currently available for com- 
parison. More data is expected in the near future for 
additional comparison. 

Figure 3 show's some typical comparison of the wake 
LIDAR data with the present calculations. These sim- 


ulations were performed using 61x101x101 grid points 
and using the temporal method. The initial condi- 
tions and computation domain was chosen to match 
the data for B737. For comparison, identical location 
(i.e., plume age), domain size and contour distribution 
is used in these figures. The chosen location is approx- 
imately Ikm downstream of the engine exhaust (at a 
plume age of 5 seconds). Some interesting features can 
be discerned from this comparison. Note that, the nu- 
merical study enforced symmetry along the aircraft 
center-plane primarily to reduce the computational 
cost (see Figure la). As a result, the numerical result 
shows a symmetric image (obtained by mirror imaging 
the data) whereas the experimental data shows some 
asymmetry. In spite of this difference there is quite 
a bit of similarity in the predicted and the measured 
vortex structure. These features are discussed below. 

First, the scale (size) and the location of the com- 
puted structure is in good agreement with data. The 
separation distance between the vortex pair at this 
demonstration location is predicted quite accurately 
even with symmetry boundary condition. Second, the 
structure of the vortex with its trailing tail (which is 
the jet plume as it gets entrained) is quite similar. In 
both the experimental data and in the present simu- 
lation, the concentration variation in the core of the 
vortex is very small. However, there is relatively a 
large variation in the outer region of the vortex. This 
implies that the core of the vortex has a stabilizing 
effect and inhibits the dilution of the passive tracer. 
The turbulent small-scale structures seen in the exper- 
iments are captured by the present study. On the other 
hand, the earlier comparison reported using the I NI- 
WAKE method 16 did not yield this type of detail on 
the vortex structure and showed a scalar field that was 
more concentrated and coherent than the experimental 
data since the small scale mixing effects are neglected. 
The more reasonable agreement between the present 
LES and the data suggest a potential advantage of 
carrying out full 3D LES where the unsteady mixing 
effects are explicitly simulated. 

We expect to carry out additional comparison with 
LIDAR data when it becomes available in the near 
future. 

4.2 Plume- Vortex Interaction behind a 
Boeing 747 

The LES model was also used to study the B747 
wake. The key differences between the B737 and B747 
are in the appropriate scales, flight conditions and the 
presence of two engine plumes for the B747. These 
simulations were carried out using both temporal and 
spatial methods to facilitate direct comparison. The 
temporal simulation captures a region of plume-vortex 
evolving in time which corresponds spatially (via a ve- 
locity transformation) to a location that moves (with 
increasing time) in the downstream direction. A grid 
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resolution of 121 x 97 x 97 is used for both spatial and 
temporal simulations. For some simulations (of an ex- 
tended domain) a resolution of 201 x ill x 111 was 
employed. The domain size (unless otherwise noted) 
is 24s in streamwise direction and 2.5s in other two 
directions. Here, s is the aircraft's semispan. Only 
one-half of the wake regime is simulated by employing 
the symmetry boundary conditions at the wake center 
plane. 

4.3.1 Comparison with experimental data 

Figure 4 shows a comparison between present spatial 
and temporal simulations of the passive tracer con- 
centration on a cross-section of the plume/wake at a 
distance of 30 semi-spans downstream of the aircraft. 
The plume age is about 5 seconds. As seen earlier in 
Figure 3. large-scale structures are seen in both the 
data and the LES. The wing-tip vortex and entrain- 
ment tail of the jet plumes can be clearly seen in the 
present computation. Again, the concentration varia- 
tion in the core of the vortex is very small even for the 
twin plume interaction for the B747. More, turbulent 
small-scale structures are seen for the B747 due to the 
interaction between twin plumes with the vortex. 

A key difference is seen between the temporal and 
spatial data in the structure of the plumes as it gets en- 
trained into the vortex. It appears that a large portion 
of the inboard plume is still very distinctly separated 
from the vortex for the spatial case whereas in the tem- 
poral case the inboard plume is also getting entrained 
into the vortex. The implication of this difference is 
discussed in the next section. 

The descending vortex velocity is relatively steady 
within the computational range. The computed de- 
scending vortex velocity for spatial simulation and 
temporal simulation are 0.90 m/s and 0.86 m/s, re- 
spectively. From the data by Teske et a/., lb the de- 
scending vortex velocity is around 0.93 m/s between 
about 15 and 40 semi-spans which is in agreement with 
present simulations. 

Figure 5 shows the comparison of the averaged SO 2 
scalar concentration decay with increasing plume age. 
Here, SO 2 is a marker species present initially only in 
the plume and therefore represents the transport and 
entrainment of the plume into the vortex. The current 
result has been normalized by the SO 2 concentration 
value in Teske et at. 16 at plume age 1 second in order to 
compare the dilution process between the two models. 
This was required since the present simulation starts 
with a different initial scalar concentrations when com- 
pared with Teske et al. 16 The reason is that there is 
no experimental data at a plume age of 1 second and 
the inflow setup in Teske et al 16 and the present LES 
are not similar. The averaged passive scaler concen- 
tration is defined as the volume average of the cells 
which has at least one percent of the local maximum 
concentration. This definition was also used by Teske 


et al. 16 

The decrease in SO 2 concentration is mainly due 
to mixing with ambient air because the reaction rate 
for SOo destruction is small. The current spatial and 
temporal simulations are in good agreement, but there 
is a discrepancy between present study and t he re- 
sult obtained using UNIWAKE. lh The current result 
shows a much larger initial decay rate than the UNI- 
WAKE result. The experimental measurements by 
Arnold et al. 32 also showed that a lower value of SO 2 
at a distance of 2km downstream of the plane than the 
prediction by Teske et a/. 16 The current prediction in 
that region shows a better agreement with data. As 
Teske et a/. 16 noted, the limitation of the turbulence 
model used in UNIWAKE may be contributing to this 
discrepancy. In the present LES with the dynamic 
sub-grid model the turbulent viscosity is dynamiclv 
computed and adjusts to varying strain rate of the 
mean flow (i.e., different age of the plume). Thus, 
better comparison with the experimental data is ex- 
pected and can be observed in this figure. It can also 
be argued that the more accurate resolution of the 
turbulent small-scale structures in the LES improves 
the prediction since turbulent diffusion plays a large 
role in scalar mixing. In contrast, UNI WAKE predic- 
tions ignore the small scale mixing effects. This may 
contribute to the under-prediction of average scalar 
concentration decay. 

4.3.3 Detramment of the exhaust plume. 

Figures 6a, 6b and 6c show respectively, three in- 
stantaneous vorticity magnitude contour plots in the 
x-z plane at y=0 from the spatial simulation. The two 
jet plumes develop almost independently before the 
first 4 semi-spans. After 4 semi-spans, the outboard 
jet plume starts to get entrained into the wing-tip vor- 
tex, and the jet plume is deflected towards the wing-tip 
vortex. Further downstream, the two jet plumes break 
up into smaller structures. This is characterized by 
patches of concentrated vorticity which are surrounded 
by flow with much lower vorticity. 

The high level of intermittency seen in the flow field 
makes Reynolds averaged approach inapplicable since 
it cannot capture this effect. It can be seen that the 
broken portions of the outboard jet plume are de- 
flected and entrained into the wing-tip vortex. On the 
contrary, due to the strong turbulent motion and rela- 
tively weak influence of the wing-tip vortex, a portion 
of the inboard jet plume gets detrained from the wing- 
tip vortex. Observation of contrails from widebody 
planes such as B747 show that there are two contrails 
behind the aircraft, one that is the normal vortex core 
contrail and another that is distinctly separated from 
it. As Gerz and Eliret 15 pointed out. irregular roll- 
up of wing-tip vortex, turbulent mot ion and buoyancy 
may result in detrainment. Gerz and Eliret u failed 
to capture detrainment in the near field by including 
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buoyancy effect. Also, the detrainnient captured by 
Quackenbush at e/. 31 was only in the far field. Present 
simulation suggests another reason for detrainnient. 
This is the effect of turbulent mixing which breaks up 
the inboard jet plume before it gets entrained into the 
vortex. This in turn forms the detrained contrail as 
seen in Figures 6a-c. 

The instantaneous NO mass fraction field corre- 
sponding to Figures 6 are shown in figures 7a, 7b and 
7c, respectively. In these simulation, the species do 
not to have a strong effect on the fluid flow because the 
reactions release almost no heat. As a result, the corre- 
lation between the NO mass fraction and the vorticity 
magnitude is very strong. Since the reaction rate for 
NO species is relatively small, the NO mass fraction 
distribution can be used to represent the overall plume 
mixing process. The entrainment and detrainnient of 
the NO species are thus similar to the vorticity field. 
It is worthwhile to point out that previous numerical 
studies 16 ’ 31 failed to capture the plume breakup and 
instead showed a smoothly varying plume species field. 
This can result in error in entrainment prediction. 

Figure 8 shows a 3-D snapshot of the instantaneous 
vorticity magnitude isosurface with NO mass fraction 
contour. The sequence of the plume-vortex interaction 
starts with the break-up of the jet plumes. By the mid- 
dle point (x=12s), the two jet plumes have broken up 
completely and part of the outboard jet plume starts 
to get entrained into the vortex. At the end of the 
interaction, a major part of jet plumes has been en- 
trained into the wing-tip vortex. However, a portion 
of the inboard jet plume is detrained from the vortex. 

The simulations by Gamier et a/. 12 showed that 
the outboard jet plume scalar concentration decays 
faster than the inboard jet plume. They attributed 
this phenomenon to the wing-tip vortex effect. How- 
ever, the current simulations show that both inboard 
and outboard plumes have almost the same decay rate 
even though the outboard plume is deflected toward 
the wing-tip vortex. When the outboard plume in- 
teracts with the wing-tip vortex, due to the shear in 
the flow, strong turbulence is produced. As a result, 
scalar mixing is enhanced and therefore, faster decay 
may be expected. However, when part of the plume is 
entrained into the core of the vortex, due to the inher- 
ent stability of the vortex core region, scalar mixing is 
slowed down. These two contradicting effects result in 
similar mixing rate for the outboard and the inboard 
plumes. 

Another explanation for the results obtained by 
Gamier et al. 12 is related to the well-known over pre- 
diction by the k — e turbulence model. The over 
prediction of the overall S0 2 decay rate by Gamier 
et al. 12 is also attributed to k - e turbulence model. 

Figures 9a and 9b show respectively, two instanta- 
neous cross-section flow field of the vorticity magni- 
tude contour plots at the downstream location 12s and 


24s (s is the semi-span) from the temporal simulation. 
Figures 10a and 10b show the corresponding vorticity 
contour from the spatial case. There is great variety 
in the instantaneous cross section but these figures are 
very representative. At 6s. the two jet plumes become 
turbulent and the interaction between the outboard 
jet plume and wing-tip vortex starts (not shown). At 
12s, the two jet plumes begin interact with each other 
and merge. However, some part of the outboard jet 
plume gets entrained into the wing-tip vortex. For the 
spatial case, at 18s (not shown), a large portion of 
the outboard jet plume is entrained into the wing-tip 
vortex while a major portion of the inboard jet plume 
stays separated. 

The detrainnient of the inboard plume continues as 
shown at x=24s and is more apparent in the spatial 
case than in the temporal case. Although there exist 
highly turbulent structures around the wing-tip vor- 
tex, the core of wing-tip vortex remains coherent and 
laminar. This overall observation of the plume disper- 
sion is quite similar to the results obtained by Sykes 
and Henn 33 and Sykes et al., 34 except for the addition 
of wing-tip vortex in present study. 

Comparison of the temporal (Figs. 9) and the spa- 
tial (Figs. 10) cases show some significant differences 
even through the overall entrainment process is qual- 
itatively similar. For example, in the temporal cases, 
the entrainment process is not very well defined in 
comparison to the spatial case. The detrained plume 
for the spatial case is also much more coherent than 
in the temporal case (also shown in figure 4). The 
inability of the temporal simulation to capture the de- 
trainment process is due to the fact that detrainnient 
occurs right behind the aircraft where the periodic as- 
sumption is not satisfied. This observation suggests 
that to capture the observed detrainment of the plume, 
spatial simulations are required. 

4.2.3 Comparison between spatial and temporal 
simulations 

Some differences between temporal and spatial sim- 
ulations have been noted above. Namely, the spatial 
simulation is able to capture detrainment effects while 
the temporal case could not. This was attributed to 
the fact that in the early stages of plume- vortex inter- 
action periodic assumption in the streamwise direction 
may not be appropriate. Here, to further quantify 
the results additional entrainment quantities are com- 
puted and compared. Figure 11 shows the detrainment 
ratio from the vortex core. The detrainment ratio is 
defined as the percentage of representative species NO 
outside the wing-tip vortex core which has a radius 
about 4m. The vortex center is located by finding 
the point where the minimum axial vorticity occurs. 
Initially, the detrainment rate is 1.0. The spatial and 
temporal simulations show similar trends however, the 
spatial simulation data has some oscillations. The os- 
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dilations of the spatial data could be attributed to 
the random wandering of wing-tip vortex which has 
been observed in the experiment by Devenport et a/. 35 
Thus, it would appear that spatial simulation is able 
to capture this random wandering effect. 

As the plume begins to get entrained into the vortex 
core, the detrainment ratio decreases. This decrease is 
large in the first 12s and then levels off. This indicates 
that part of the plume has overcome the attraction of 
the wing-tip vortex after 12s. This also can be seen 
in figure 7a, 7b and 7c where there is an apparent 
separation of jet plumes after 12s. The detrainment 
reaches a final value of around 0.87. Therefore, about 
13 percent of jet plume is entrained into the core of 
the wing-tip vortex. 

In the vortex core, the temperature is about 8 de- 
gree higher than the ambient temperature. Note that 
without entrainment of the hot jet plumes, the temper- 
ature inside the wing-tip core should be a few degree 
lower than the ambient temperature due to the lower 
pressure inside the vortex core. Due to entrainment, 
the temperature of vortex core is increased. The sev- 
eral degrees' increase of temperature is critical when 
chemical reactions are included. This difference in 
temperature indicates that the entrained plume un- 
dergoes a different chemical reaction path, e.g., via the 
heterogeneous reactions when compared to the portion 
that is detrained from the vortex. The previous stud- 
ies 16,31 showed only a fully entrained plume without 
any details of the structure inside the vortex. 

It is interesting to track the different dilution pro- 
cesses of plume species inside and outside the vortex 
core since there is a lower pressure zone inside the vor- 
tex core and the fluid field is more stable inside the 
vortex core. It is well known that the conditions for 
contrail formation are very subtle and a little differ- 
ence in the thermodynamics and fluid dynamics can 
result in widely different contrail formation. Figure 12 
shows the maximum NO mass fraction inside and out- 
side the vortex core versus downstream distance. Both 
spatial and temporal results are compared. Inside the 
vortex core, the maximum NO mass fraction is almost 
constant after it reaches the maximum. While, outside 
the vortex core the maximum NO mass fraction drops 
very fast in the first 12s. The reason is that turbulent 
diffusion is very large outside the vortex core due to the 
enhanced mixing effect due to strong turbulence where 
jet plume breakdown occurs. On the other hand, the 
vortex core is relatively stable due to the swirling effect 
of the wing-tip vortex. After 18s, when the turbulence 
begins to die down due to dissipation, the NO concen- 
tration begins to level off. Comparisons of the spatial 
and temporal results show that on the outside max- 
imum NO mass fraction agrees very well. However, 
inside the vortex core, the temporal simulation shows 
a higher peak than the spatial simulation (by around 
257c) . 


Figure 13 shows the axial variation of NO mass 
fraction at vortex center. Although the spatial and 
temporal simulations show similar trend, the spatial 
simulation shows a much higher mass fraction than 
the temporal simulation. From figure 12 and 13, one 
can conclude that the species distribution inside the 
vortex core predicted by the temporal simulation is 
more spread out around the core center than the pre- 
diction by the spatial simulation. This implies that 
the reaction rate in the temporal simulation is slower 
than that in the spatial simulation. 

Figure 14 shows the variation of SO 3 mass fraction 
at vortex center with downstream distance for both 
spatial and temporal simulations. In the spatial simu- 
lation, SO3 mass fraction first increases and then levels 
off near 24s. On the other hand, in the temporal sim- 
ulation, SO 3 mass fraction reaches a peak at around 
12s and then decreases. There are significant differ- 
ence in the magnitude of peak SO 3 . The spatial data 
shows 10” 8 where as the temporal data shows 10“ u . 
The wide discrepancy between temporal and spatial 
results has important implications for heterogeneous 
reactions and for H 2 SO 4 aerosol formulation. This is 
briefly discussed below. 

Note that SO 3 is produced first by the oxidization 
of SO 0 and then primarily removed by heterogeneous 
reaction on condensed soot particles. The condensa- 
tion inside the core may be small due to the lower 
pressure. As a result, the removal rate of SO 3 is small 
so that SO 3 mass fraction can keep increasing down- 
stream. The magnitude of the SO 3 mass fraction is 
thus determined by the rate of SO 2 oxidization reac- 
tion which is strongly affected by SO 2 mass fraction. 
As figures 12 and 13 suggest, the spatial data has a 
higher species concentration around the vortex center 
than the temporal data. This explains the higher SO 3 
concentration in the spatial simulation. These differ- 
ences between the spatial and temporal simulations 
show the necessity for spatial simulation. 

The implication of the higher concentration of SO 3 
is that it could result in an increase in the H 0 SO 4 
aerosol formation. Although this feature has not been 
included here, the present spatial LES prediction ap- 
pears to partly explain the discrepancy between ex- 
perimental data and earlier numerical studies. As 
noted earlier, 36 data suggests a 45 times higher H 0 SO 4 
aerosols concentration in the vortex- plume than the 
earlier calculations (using SPF-II type of code). The 
present prediction of higher SO 3 concentration ob- 
tained by doing a spatial simulation (in contrast to 
temporal or parabolic space marching) seems to sug- 
gest that this discrepancy may be partly due to the 
numerical approach. However, to confirm this we need 
to carry out these simulation including H 0 SO 4 binary 
nucleation and coagulation aerosol models. This effort 
is underway and will be reported in the near future. 

Some integrated quantities are also analyzed as fol- 
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lows. The mass of entrained NO into the wing-tip 
vortex is shown in figure 15. Temporal simulation 
shows a larger amount of NO is entrained than the 
spatial simulation. Since the relative reaction rate for 
NO is small, .VO entrainment can be used to repre- 
sent the entrainment process. The difference between 
the spatial and temporal simulations can be due to the 
axial velocity deficit of the jet plume. The deviation 
between the two types of simulations begins around 
x = 12s, where most of outboard engine plume is en- 
trained into wing-tip vortex and the inboard engine 
plume splits into two parts. One part eventually gets 
entrained into the wing-tip vortex while the second 
part becomes detrained from the wing-tip vortex. The 
chemical processes that result due to the A O concen- 
tration will be different in these two regions. 

Other integrated quantities of species also show sig- 
nificant deviation when the two simulations are com- 
pared. Figure 16 shows the mass entrainment of b 03 
into the wing-tip vortex. The result of spatial simula- 
tion shows a higher peak than the temporal simulation 
and shows an overall larger entrainment. This leads 
to the conclusion that the temporal simulation is not 
suitable for the early plume- vortex interaction region 
because there exists no periodicity in the interaction 
region. The discrepancy in the SO3 concentration in 
the vortex is particularly bothersome as noted ear- 
lier. The increased number of aerosol measured in the 
wake 36 and the potential increase in the condensed 
sulfuric acid (due to the reaction between SO3 and 
condensed water) forming on these particles have been 
identified as one process that could impact long term 
global atmosphere chemical balance. 

5 Conclusions 

In this study, fluid mixing and chemical reactions 
in the near-field wake of an aircraft has been inves- 
tigated. The plume-vortex regime is studied using a 
3D parallel LES code. Simulations of engine plume 
interaction with the wing-tip vortex of B737 and B747 
have been carried out and some comparison with data 
have been used to demonstrate the ability of the LES 
method. Inclusion of detailed gas phase kinetics and 
simple heterogeneous processes have been carried out 
to evaluate the effect of turbulent mixing on the chem- 
ical processes in the plume-vortex regime. Finally, to 
determine the impact of temporal assumptions (often 
used in such flow simulations) comparison with spatial 
simulations were carried out. 

Many features of the computed flow field show rea- 
sonable agreement with experimental data. The en- 
trainment and detrainment of engine plumes into the 
wing-tip vortex were numerically captured in the spa- 
tial simulation and analyzed. The impact of fluid 
mechanics on the chemical process has been estimated. 
Results indicate that a significant difference between 
spatial and temporal simulation exists (especially in 


the predicted SO 3 concentration). This has important 
implication for the prediction of sulfuric acid aerosols 
in the w'ake and may partly explain the discrepancy 
between past numerical studies and the data. Anal- 
ysis of the results also show that spatial simulation 
is more suitable for the near field interaction process. 
For example, it was determined that detrainment (ob- 
served in the B747 contrail) was only captured using 
the spatial model. 
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Fig. 2 Comparison of averaged S0 2 mole fraction 
decay for different turbulence inflow conditions. 
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Fig. 4 Comparison of scalar concentration be- 
tween spatial and temporal simulations 
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Fig. 9 Vorticity magnitude contours at different 
locations in the temporal simulation. 



Fig. 13 NO mass fraction at the vortex core center. 
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Fig. 10 Vorticity magnitude contours at different 
locations in the spatial simulation. 



Fig. 11 Variation of the detrainment ratio of NO 
species with the downstream distance for B747. 





Fig. 12 Maximum NO mass fraction inside and 
outside the vortex core (vortex core has a radius of 
4m). 



Fig. 16 SO3, mass entrainment into the vortex core. 
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Abstract. Species mixing and chemical reactions in the near-field wake ot an 
aircraft have been investigated using a large eddy simulation (LES) model. Simulations 
of the interaction between engine plumes and the wing-tip vortices of Bi-b and 
B747 have been carried out and comparisons with available data have been used to 
demonstrate the ability of the LES method. The effect of turbulent mixing on the 
chemical processes in the plume-vortex regime has been studied using detailed gas 
phase and simple heterogeneous kinetics. Finally, to determine the impact ot temporal 
assumptions (often used in such flow simulations) comparisons with spatial simulations 
were carried out. 

Many features of the computed flow field show reasonable agreement with 
experimental data. Analysis of the results also shows that spatial simulation is more 
suitable for the near field interaction process. The entrainment of engine plumes into 
the wing-tip vortex were numerically captured in the spatial simulation and analyzed. 
However, detrainment (observed in the B747 contrail) was only captured using the 
spatial model. Results indicate that a significant difference between spatial and tempora 
simulation exists, especially in the predicted SO 3 concentration. This has important 
implication for the prediction of sulfuric acid aerosols in the wake and may partly 
explain the discrepancy between past numerical studies and measured data. 



1. Introduction 

The wing-tip vortices shed from an aircraft s wings are known to pose a clanger 
for aircrafts in the wake and are considered a major hazard near airports. Another 
consequence that is perhaps less appreciated is their ability to entrain the hot exhaust ot 
the engine during cruise flight in the upper atmosphere. The resulting species-rich and 
relatively long-lived contrails have been of increasing concern and a focus of research 
in recent years because the species entrained into the vortex wake can react with the 
ambient species, thereby, changing the local chemical balance. These studies have been 
motivated bv the projected increase in future subsonic flights in the upper troposphere 
by the proposed Advanced Subsonic Transports (AST) operating in the transatlantic 
and Pacific rim flight corridors and by the possibility of future supersonic flights in the 
lower stratosphere by the High Speed Civilian Transport (HSCT). 

The near-field region of the aircraft's wake can be divided roughly into two parts: a 
plume regime and a plume-cortex regime. The plume regime is identified as the initial 
region behind the aircraft where the engine exhaust plume grows and mixes with the 
ambient air. This regime is typically very short, approximately 1 sec after exhaust 
or around 200 m downstream of the engine exit plane. Further downstream is the 
plume-vortex regime where the exhaust plume is entrained into the rolled-up wing tip 
vortices. This regime is around 1 - 2 km long and approximately 10 sec downstream 
of the exhaust. Due to the relatively high temperature of the engine exhaust and the 
species content, significant reactions can take place in the plume-vortex regime. The 
region downstream of the plume- vortex region is the wake region where the vortex wake 
breaks up due to the Crow instability and mixes with the ambient air. This region is 
very large extending over 20 km and the mixing process can take many days depending 
upon the ambient conditions. 

The near-field region, especially the plume- vortex regime, is of interest here. In 
particular, we are interested in the chemical processes that occur in this region. Key 
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species of interest are nitrogen oxides {NO x ), carbon monoxide (CO), methane (CH A ). 
non-methane hydrocarbons ( HC ) and soot particles. These trace species pla\ an 
important role in the atmospheric photochemical and radiative processes. For example, 
the concentration of tropospheric ozone, an important greenhouse gas. is dependent 
on transport across tropopause and by chemical production and destruction in\ol\ing 
reactions of .VO r , non-methane hydrocarbons, carbon monoxide, and hydrogen oxides 
( HO x ). Furthermore, the radiative balance of the atmosphere can be affected by 
the accumulation of carbon dioxide and water and by the aircraft-induced changes 
in aerosols and clouds. In addition, the combustion aerosols (soot) and the binary 
H 2 S0 a /H 2 0 droplets in the plume can be efficient cloud condensation nuclei and serve 
as sites for heterogeneous reactions. 

A varietv of models have been developed to assess the effect of aircraft flights on the 
atmosphere. These models generally fall into two categories: small-scale and global 

models. "Small-scale” models treat the near-field regime of the exhaust wake by ignoring 
processes that operate on larger time and space scales. Global models [Schumann. 1996]. 
on the other hand, attempt to represent the large spatial scales and ignore, or represent 
very crudely, the small scales processes. To obtain some sort of continuity, some of 
the global models [Schumann. 1996] use the engine exhaust conditions or the data at 
the end of the plume-vortex regime as input to their large-scale models. However, as 
shown recently [Brown et ai. 1996; Wang and Chen. 1997; Menon and Wu. 1998], the 
concentration of species (e.g., chemical radicals OH. H0 2 and soot) can undergo very 
rapid changes in the near- field. This can significantly modify the input conditions for 
the far field simulations. Thus, to obtain accurate prediction of the global chemical 
balance due to aircraft emission, the chemical processes in the plume-vortex regime are 
very important. 

The hydrodynamics of the plume- vortex regime is also quite complex. For example, 
the competition between the vortex centripetal force and gravitation buoyancy of 



hot plume can lead to a different final altitude of deposition of emissions [Shake -Lge 
et al. . 1993]. This implies that evolution of the plume occurs in a highly unsteady 
manner. However, many of the models [Brown et al., 1996: Wang and Chen. 1997] are 
steady-state models and therefore, are unable to account for the unsteady mixing effects. 

To study unsteady jet-vortex interactions, unsteady methods have also been used. 
For example, direct numerical simulation [Risso et al., 1996] and large eddy simulation 
(LES) of trailing vortices [ Ragab and Sveedhar , 1994a] showed that there exists strong 
3-D instability in vortex wakes with an axial deficit. The vortex instability results in 
vortex-stretching which may have a substantial effect on the mixing of jet plume with 
the ambient air. Ger: and Ehret [1996] used LES to simulate the vortex roll-up process 
of the wake of a B-747 one second after exit. They showed that temperature, relative 
humidity peak at the center of the wingtip vortex: the wing boundary layer turbulence 
and the ambient turbulence all can affect the plume- vortex interaction process. 

Unsteady simulations are computationally expensive and therefore, most past 
studies employed some simplifying assumptions. For example, only limited chemical 
effects have been investigated so far [Teske et al.. 1997] and all simulations reported so 
far employed periodic boundary conditions in the streamwise directions (i.e. simulated 
temporal evolution instead of the actual spatial evolution). Temporal approximation 
is truly justified only far downstream where the interaction between the jet plume and 
wingtip vortex has been completed [Lewellen et al., 1998]. It has not yet been shown 
what features of the flow and the accompanying chemical processes are different and/or 
modified by this temporal assumption in the near field. This issue is important since 
it has been shown that the entrainment predicted by temporal and spatial simulations 
differ markedly due to the asymmetry between the high speed and the low speed sides 
of the shear layer [Dimotakis, 1986]. This asymmetry cannot be captured in a temporal 
simulation. The impact of this asymmetric entrainment on the chemical processes in 
the plume- vortex interaction regime has not been addressed so far. 
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This paper will address many of these issues and will attempt to quantify the 
impact of simulation methodology (i.e., spatial or temporal) on the prediction ot the 
species distribution. Furthermore, the effect of the fluid dynamic interaction process on 
the chemical kinetics and on the concentration of species such as A'0, 503 and 0$ in 
the plume- vortex regime will be addressed in the presence of both gas and heterogeneous 
kinetics. 


2. Simulation Model 


The formulation of the simulation model is briefly summarized below. More details 
are given elsewhere [Menon et al . , 1996; Kim and Mtnon. 1995; Kelson and Me non. 
1998] and. therefore, avoided here for brevity. 

The fluid dynamic LES equations are obtained by spatially filtering the compressible 
Xavier Stokes equations. The resulting equations contain unknown subgrid terms such 
as the stress in the momentum equations. T- gs = p[upij ~ the enthalpy flux and 

the viscous terms in the energy equation. H t S9S = p[Hu, — Hu,\. a- gs - [ iijT ,, — 
and the scalar correlations. h\ 9$ = TY k - TYk in the equation of state. Here, tilde 
and bar indicate filtered variables, tq, p, f tJ and H are. respectively, the resolved 
velocity components, the density, the viscous stresses and the total enthalpy per unit 
volume. To close these subgrid terms, a model equation for the subgrid kinetic energy 


k 3gs — u,u, — u,iii : 


dpi ; 3gJ 

dt 


+ l-(pu i k t n 

axi 


d , _ dk sg3 , 

d7, (pv '- 

-r-fditj/dxj 


-c,p(k 393 f /2 /K 


( 1 ) 


is also solved along with the LES equations. Here, v t = Cvsfk 3 * 3 A is the subgrid eddy 
viscosity and A is the grid size. The coefficients c v and c c are computed locally (in space 
and time) during the simulation. With k 3gs known, the subgrid terms are approximated 



as T*f* « -2 pv t (Sij - lSkkSij) + | pW6 u , Hr = and <7**' = urff. 

Here. Pi' t % 1 is the turbulent Prandtl number and is the resolved rate-of-stiain 
tensor. S tJ = |(H*- + |^). All third-order subgrid correlations are neglected in this 
closure for simplicity. 

Closure of the species equation requires closure of the subgrid species transport 
term p\\ uj — pVl-Uj ■ A subgrid diffusivity model (as the eddy viscosity model given 
above) is not considered an excellent choice since for combustion to occur, species must 
first mix at the small scales and then molecularly diffuse. These processes occur at the 
unresolved scales and therefore, a global eddy diffusivity model for the subgrid species 
flux that ignores the subgrid processes, requires some justifications and can be used only 
under certain conditions. In the present study, the mixing process is dominated by the 
large-scale transport of the plume species into the wing-tip vortices. Furthermore, the 
time-scales for the subsequent chemical reactions are quite large (i.e.. slow chemistry). 
Thus, the fast time small-scale mixing effects can be considered completed before the 
chemical reactions occur. For these reasons, an eddy diffusivity model for the species 
subgrid flux is considered a reasonable first-order approximation. With these comments 
in mind, the subgrid species flux is approximated as: 

— r= — -v- - _ P Vt dYk r» 

P*k u j P*kUj — q c T Q r . 

where Sc T is the turbulent Schmidt number set to unity for the present study. 


3. Numerical Method 

An explicit finite-volume scheme that is fourth order accurate in space and second 
order accurate in time was used in the present study. The details of the numerical scheme 
has been reported elsewhere [Nelson, 1997; Nelson and Menon , 199S] and therefore, 
avoided here for brevity. The current approach is similar to the implementation 
described by Garnet and Estivalezes [1995] except that the fluxes are modified so that a 
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truly fourth-order accuracy in space is obtained at least on an uniform grid. 

The present code has been parallelized using the standard Message Passing Interlace 
(MPI) library. The computational domain is decomposed into subdomains and assigned 
to different processors which have an identical copy of the numerical algorithm. For the 
current scheme, a five point stencil is used and therefore, two rows of variables must be 
exchanged in either direction. Timings analysis of this code on various parallel machines 
(e.g., Cray T3E. SGI Origin 2000) show that the code achieves nearly linear speedup as 
the number of processors is doubled. However, depending upon the total number of grid 
points employed there is an optimum number of processors beyond which the message 
passing overhead overtakes the computing time in each processor s domain. 

For a 121 x 97 x 97 computational grid a typical simulation (using lb reactive 
species) on the Cray T3E takes around 5000 single processor hours to obtain data 
sufficient for statistical analysis (-5 flow through time). The CPU time per time step per 
grid point was around 3E-6 sec on 12S processor T3E (or 3.6 sec per time step). 

Two types of problems are studied in the present study. Figures la and 11) show, 
respectively, the plume-vortex domain behind a Boeing 737 and a Boeing 717 aircraft. 
The primarv differences are in the span, the number of engine plumes and the exhaust 
conditions. These two configurations were chosen primarily because there is some data 
available for comparison. 

For slip or symmetry conditions, values are prescribed for ghost-cells such that 
gradients normal to the face at the bounding face are zero. Periodic boundary conditions 
are prescribed in the streamwise direction for the temporal simulations. 

For the spatial cases, the inflow and outflow boundaries are more complicated 
and must be carefully implemented. We employ a method that attempts to provide 
reasonably accurate inflow turbulence statistics [Nelson. 1997: Pomsot and Leif. 1992]. 
The inflow is critical for the spatial case since it determines the turbulence entering 
the computational domain. A turbulent inflow field is generated based on a prescribed 
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energy spectrum and turbulence intensity and then superimposed on the mean held. 
The turbulence spectrum chosen is of the form E(k) = C'k 4 t.vp{ —2{k/k p ) 2 ) used by Let 
tt al. [1991]. where constant C is chosen as turbulence intensity level. The turbulence is 
divergence free and non-periodic in time. The field is stored as 3D array from \\ hich a 
2-D plane is introduced at the inflow at every time step. At least 10 different arrays of 
97 x 97 x 97 are stored and to ensure that the inflow field is not correlated, at randomly 
chosen instants, the phase in the Fourier modes is randomized. 

The outflow boundary is non-reflecting characteristic boundary conditions [Poinsot 
and Ltle. 1992]. The amplitudes of the outgoing waves are computed from the 
interior points and the amplitude of the incoming wave is computed by specifying the 
background pressure. Details of this type of outflow condition are given in [Pom.-ot and 
Ltle. 1992], 

The wing-tip vortex is modeled by the Lamb-Oseen vortex [.S off in an. 1973] with 
an axial velocity deficit. This model has been extensively used [Ragab and Srttdbar. 
1994a. b: Quacktnbush tt al.. 1993. 1996; Risso tt al.. 1996: Gamier tt al. 1996] for 
earlier studies. The location, core radius, maximum circulation data of the wing-tip 
vortex and the engine location and exit conditions (Table 1 and 2 summarize these 
conditions) are chosen according to the aircraft configurations. 

The gas phase chemical kinetics studied in the present paper are identical to the 
chemical mechanism in [Brown tt al.. 1996] except that in the present case, we employ a 
reduced reaction mechanism that shows an excellent agreement with the full mechanism 
[ Wang and Chen. 1997; Mtnon and Wu. 1998]. A total of 15 species were simulated in 
the present study. The heterogeneous kinetics studied are very simplistic and consist of 
condensation of nitric and sulfuric acid on the soot particles. The details of this model 
are given in [I Vang and Chen. 1997; Mtnon and Hu, 1998] and avoided here for brevity. 


Table 1. 1 
Table 2J 
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4. Results and Discussion 

In this section, we summarize the results obtained in the present study. There 
are three primary objectives of this study. First, comparison with available data is 
carried out to demonstrate the ability of the simulation model to accurately capture the 
evolution of the plume-vortex regime. Second, the results obtained using temporal and 
spatial methods are compared to identify what features are not predicted accurately 
using temporal simulations (as noted earlier, nearly all previously reported 3D studies 
employed periodic conditions in the streamwise direction). The final objective is 
to determine the impact of turbulence transport on the formation of SO r and the 
destruction of ozone. 

Before discussing the details of the plume-vortex interaction, the effect ot different 
inflow turbulence on the flow field is briefly summarized. Figure 2 shows the effect ot 
inflow turbulence on the scalar (here S0 2 concentration is used as a scalar marker) 
decay with increasing plume age in a B747 plume-vortex region. Three different inflow 
turbulence conditions are used. Two conditions employed the same energy spectrum 
but with a 5 % and 10% turbulent intensity, respectively, while the third case used the 
same initial isotropic spectrum as the other two but was introduced into the flow field 
after evolving for a certain length of time with a final turbulent intensity ol • 7%. Thus, 
in the third case, the inflow evolved to a more realistic turbulent isotropic state whereas 
in the first two cases, the field is isotropic but is not realistic turbulence. 

It can be seen that the cases with the same turbulent intensity yield almost identical 
results. This means that the effect of evolving of the inflow turbulent field is small, at 
least for the global measure used here. (Note that such global measures are of interest 
here). However, it can be seen that turbulent intensity does have an effect on the flow 
field. This implies that the inflow turbulent intensity must be carefully chosen. The 
turbulence intensity in the vortex and jet plume varies with flight conditions such as 
climate, altitude and ambient flow field. A turbulent intensity between o% and 10% is 
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very common in the typical flight path of interest [Ragab and Sreedhar , 1994b]. Thus, 
simulations with 5% and 10% inflow turbulence may be used to quantify an upper and 
lower bound for the species distributions in actual flights. 

Interestingly, the descending vortex velocity computed from these three cases 
showed similar results. However, this is understandable because the descending vortex 
velocity is mainly determined by the overall vortex circulation and small scale turbulence 
exerts negligible effects on the circulation. The numerical simulation of q- vortex by 
Ragab and Sreedhar [1994b] also demonstrated this feature. 

4.1. The Boeing 737 wake 

Figure 3 shows some typical comparison of the wake LIDAR data [Poole e t a!.. 
1996] with the present calculations. These simulations were performed using 61x101x101 
grid points and using the temporal method. The initial conditions and computation 
domain was chosen to match the data for B737. For comparison, identical location 
( i.e. . plume age), domain size and contour distribution is used in these figures. The 
chosen location is approximately l kin downstream of the engine exhaust (at a plume 
age of 5 seconds). Note that the numerical study enforced symmetry along the aircraft 
center-plane primarily to reduce the computational cost (see Figure la). As a result, 
the numerical result shows a symmetric image (obtained by mirror imaging the data) 
whereas the experimental data shows some asymmetry. In spite of this difference there 
is quite a bit of similarity in the predicted and the measured vortex structure. These 
features are discussed below. 

First, the scale (size) and the location of the computed structure is in good 
agreement with data. The separation distance between the vortex pair is also predicted 
cjuite accurately even with symmetry boundary condition. Second, the structure of 
the vortex with its trailing tail (which is the jet plume as it gets entrained) is quite 
similar. In both the experimental data and in the present simulation, the concentration 
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variation in the core of the vortex is very small. However, there is relatively a large 
variation in the outer region of the vortex. This implies that the core of the vortex 
has a stabilizing effect and inhibits the dilution of the passive tracer. The turbulent 
small-scale structures seen in the experiments are qualitatively captured by the LES. On 
the other hand, the earlier comparison reported using the UNIWAKE method [ Teskt 
et al.. 1997] did not yield this type of detail on the vortex structure and showed a scalar 
field that was more concentrated and coherent than the experimental data since the 
small scale mixing effects were neglected. The more reasonable agreement between the 
present LES and the data suggest a potential advantage of carrying out full 3D LES. 

We expect to carry out additional comparison with the LIDAR data when it 
becomes available in the near future. 

4.2. The Boeing 747 wake 

The LES model was used extensively to study the B747 wake. The key differences 
between the B737 and B747 are in the appropriate scales, flight conditions and the 
presence of two engine plumes for the B747. These simulations were carried out using 
both temporal and spatial methods to facilitate direct comparison. The temporal 
simulation captures a region of plume-vortex evolving in time which corresponds 
spatially (via a velocity transformation) to a location that moves (with increasing time) 
in the downstream direction. A grid resolution of 121 x 97 x 97 is used for both spatial 
and temporal simulations. For some simulations (of an extended domain) a resolution of 
201 x 111 x 111 w f as employed. The calculation begins from 6s and domain size of 24s 
in the streamwise direction was simulated. In the other two directions, the domain size 
is 2.5s. Here, s is the aircraft’s semispan. Only one-half of the wake regime is simulated 
by employing the symmetry boundary conditions at the wake center plane. 

4.2.1. Plume evolution Figure 4 shows a comparison between present spatial 
and temporal simulations of the passive tracer concentration on a cross-section of the 
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plume/wake at a distance of 30 semi-spans downstream of the airciaft. The plume age 
is about 5 seconds. As seen earlier in Figure 3 (for the B737). large-scale structures 
are seen in the present LES. The wing-tip vortex and the tail of the jet plumes can 
be clearly seen in the present computation. Again, the concentration variation in the 
core of the vortex is very small even for the twin plume interaction for the B747. More 
turbulent small-scale structures are seen for the B747 due to the interaction between 
twin plumes with the vortex. 

A key difference is seen between the temporal and spatial data in the structure 
of the plumes as it gets entrained into the vortex. It appears that a large portion of 
the inboard plume is still very distinctly separated from the vortex for the spatial case, 
whereas in the temporal case the inboard plume is also getting entrained into the vortex. 
The implication of this difference is discussed in the next section. 

The descending vortex velocity is relatively steady within the computational range. 
The computed descending vortex velocity for spatial simulation and temporal simulation 
are 0.90 m/s and 0.S6 m/s, respectively. From the data by Teske et al. [199<]. the 
descending vortex velocity is around 0.93 m/s between about 15 and 40 semi-spans 
which is in agreement with present simulations. 

Figure 5 shows the comparison of the averaged S0 2 scalar concentration decay 
with increasing plume age (These simulations were carried out for large domain using 
201 x 111 x 111 grid). Here, S0 2 is a marker species present initially only in the plume 
and therefore, represents the transport and entrainment of the plume into the vortex. 
The averaged concentration is defined as the volume average of the cells which has at 
least one percent of the local maximum concentration. This definition was also used 
by Teske et al. [1997]. This was required since the present simulation begins with a 
different initial scalar concentration when compared with Teske et al. [1997], Since the 
inflow in Teske et al. [1997] and the present LES are not similar, the current result has 
been normalized by the S0 2 concentration value in Teske et al. [1997] at plume age 1 
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second in order to compare the dilution process between the two models. 

The decrease in SO 2 concentration is mainly due to mixing with ambient air 
because the reaction rate for SO 2 destruction is small. The current spatial and temporal 
simulations are in good agreement, but there is a discrepancy between present study 
and the result obtained using UNIYVAKE [Teske et al . , 1997]. The current results 
show a much larger initial decay rate than the UNIWAKE result. The experimental 
measurements by Arnold et al. [1997] also showed a lower value of SO 2 at a distance of 
2km downstream of the plane than the prediction by Teske et al. [1997]. The current 
predictions in this region show a better agreement with data. As Teske et al. [1997] 
noted, the limitation of the turbulence model used in UNIWAKE may be contributing 
to this discrepancy. In the present LES with the dynamic sub-grid model the turbulent 
viscosity is dynamicly computed and adjusted to varying strain rate of the mean flow 
(i.e.. different age of the plume). Thus, better comparison with the experimental data 
is expected and can be observed in this figure. It can also be argued that the more 
accurate resolution of the turbulent small-scale structures in the LES improves the 
prediction since turbulent diffusion plays a large role in scalar mixing. In contrast. 
UNIWAKE predictions ignore the small scale mixing effects. This may contribute to 
the under-prediction of average scalar concentration decay. 

4.2.2. Detrainment of the exhaust plume Figures 6a and 6b show, 
respectively, two instantaneous vorticity magnitude contour plots in the x-z plane at 
y=0 from the spatial simulation. The two jet plumes develop almost independently 
in the first 4 semi-spans. After 4 semi-spans, the outboard jet plume starts to get 
entrained into the wing-tip vortex, and the jet plume is deflected towards the wing-tip 
vortex. Further downstream, the two jet plumes break up into smaller structures. This 
is characterized by patches of concentrated vorticity which are surrounded by flow with 
much lower vorticity. 

The high level of intermittency seen in the flow field makes Reynolds averaged 
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approach inapplicable since it cannot capture this effect. It can be seen that the broken 
portions of the outboard jet plume are deflected and entrained into the wing-tip vortex. 
On the contrary, due to the strong turbulent motion and relatively weak influence of 
the wing-tip vortex, a portion of the inboard jet plume gets detrained from the wing-tip 
vortex. Observation of contrails from widebody planes such as the B747. show that 
there are two contrails behind the aircraft; one that is the normal vortex core contrail 
and another that is distinctly separated from it. As Gerz and Ehret [1996] pointed 
out. irregular roll-up of the wing-tip vortex, turbulent motion and buoyancy may result 
in detrainment. Ger; and Ehret [1996] failed to capture cletrainment in the near field 
by including buoyancy effect (Note that, buoyancy is also included here). Also, the 
detrainment captured by Quackenbush et al. [1996] was only in the far field. Present 
simulation suggests another reason for detrainment. This is the effect of turbulent 
mixing which breaks up the inboard jet plume before it gets entrained into the vortex. 
This, in turn, forms the detrained contrail as seen in Figures 6a-b. 

Figures 7a and 7b show, respectively, two instantaneous cross-section flow fields 
of the vorticity magnitude contour plots at the downstream location 12s and 24* (s is 
the semi-span) from the temporal simulation. The downstream distance for temporal 
simulation is estimated from the temporal time by using the Taylor's hypothesis. 
Figures Sa and Sb show the corresponding vorticity contours from the spatial case. 
There is great variety in the instantaneous data but these figures are very representative. 
At 6s. the two jet plumes become turbulent and the interaction between the outboard 
jet plume and the wing-tip vortex starts (not shown). At 12s, the two jet plumes begin 
to interact with each other and then merge. However, some part of the outboard jet 
plume gets entrained into the wing-tip vortex. For the spatial case, at 18s (not shown), 
a portion of the outboard jet plume is entrained into the wdng-tip vortex while a major 
portion of the inboard jet plume stays separated. 

The detrainment of the inboard plume continues as shown at x=24s and is more 
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apparent in the spatial case than in the temporal case. Although there exist highly 
turbulent structures around the wing-tip vortex, the core of the wing-tip vortex remains 
coherent and laminar. This overall observation of the plume dispersion is quite similar 
to the results obtained by Sykes and Henn [1992] and Sykes et al. [1992], except for the 
addition of the wing-tip vortex in the present study. 

Comparison of the temporal (Figs. 7) and the spatial (Figs. 8) cases show some 
significant differences even through the overall entrainment process is qualitatively 
similar. For example, the detrained plume for the spatial case is also much more 
coherent than in the temporal case (also shown in Fig. 4). The inability of the temporal 
simulation to clearly capture the detrainment process is due to the fact that detrainment 
occurs right behind the aircraft where the periodic assumption is not satisfied. This 
observation suggests that to capture detrainment of the plume, spatial simulations may 
be more appropriate. 

4.2.3. Spatial versus temporal entrainment effects Some differences 
between temporal and Spatial simulations have been noted above. Namely, the spatial 
simulation is able to capture detrainment effects while the temporal case could not. This 
was attributed to the fact that in the early stages of plume-vortex interaction periodic 
assumption in the streamwise direction may not be appropriate. Here, additional 
entrainment quantities are computed and compared to quantify the impact of simulation 
methodology. 

Figure 9 shows the detrainment ratio from the vortex core. The detrainment ratio 
is defined as the percentage of representative species NO outside the wing-tip vortex 
core which has a radius of about 4m. The vortex center is located by finding the point 
where the minimum axial vorticity occurs. Initially, the detrainment ratio is 1.0. The 
spatial and temporal simulations show similar trends, however, the spatial simulation 
data has some variations which may be attributed to the random wandering of the 
wing-tip vortex. This wandering has been observed in experiments by Devenport et al. 
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[1996]. Thus, it would appear that spatial simulation is able to qualitatively capture 
this random wandering effect. 

As the plume begins to get entrained into the vortex core, the detrainment ratio 
decreases. This decrease is large in the first 12$ and then begins to level off. This 
indicates that part of the plume has overcome the attraction of the wing-tip vortex after 
12s. This was also seen in Figures 6a and 6b where there is an apparent separation of jet 
plumes after 12s. The detrainment ratio reaches a final value of around 0.87. Therefore, 
about 13 percent of jet plume is entrained into the core of the wing-tip vortex. 

In the vortex core, the temperature is about 8 degree higher than the ambient 
temperature. Note that without entrainment of the hot jet plumes, the temperature 
inside the wing-tip core should be a few degrees lower than the ambient temperature 
due to the lower pressure inside the vortex core. However, due to entrainment, the 
temperature in the vortex core is increased. The several degrees' increase of temperature 
is critical when chemical reactions are included. This difference in temperature indicates 
that the entrained plume undergoes a different chemical reaction path. e.g.. via the 
heterogeneous reactions (if included) when compared to the portion that is detrained 
from the vortex. The previous studies [Quackenbush et al 1996; Teske et al 1997] 
showed only a fully entrained plume without any details of the structure inside the 
vortex. 

It is interesting to track the different dilution processes of plume species inside and 
outside the vortex core since there is a lower pressure zone inside the vortex core and 
the flow field is more stable inside the vortex core. It is well known that the conditions 
for contrail formation are very subtle and a little difference in the thermodynamics 
and fluid dynamics can result in widely different contrail formation. Figure 10 shows 
the maximum NO mass fraction inside and outside the vortex core versus downstream 
distance. Both spatial and temporal results are compared. Inside the vortex core, the 
maximum VO mass fraction is almost constant after it reaches the maximum. However. 
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outside the vortex core the maximum NO mass fraction drops very fast in the first 12s. 
The reason is that turbulent diffusion is very large outside the vortex core due to the 
enhanced mixing by the strong turbulence near jet plume breakdown. On the other 
hand, the vortex core is relatively stable due to the swirling effect of the wing-tip vortex. 
After 18s, when the turbulence begins to die down due to viscous dissipation, the .VO 
concentration begins to level off. Comparisons of the spatial and temporal results show 
that on the outside maximum NO mass fraction agrees very well. However, inside the 
vortex core, the temporal simulation shows a higher peak than the spatial simulation 
(by around 25%). 

Figure 11 shows the axial variation of NO and SO 3 mass fractions at vortex center. 
Although the spatial and the temporal simulations show similar trend, the spatial 
simulation shows a much higher mass fraction than the temporal simulation. From 
Figures 10 and 11. one can conclude that the species distribution inside the vortex core 
predicted by the temporal simulation is more spread out around the core center than the 
prediction by the spatial simulation. This implies that the reaction rate in the temporal 
simulation is slower than that in the spatial simulation. 

In both the spatial and the temporal simulations. SO 3 mass fraction first increases 
and then peaks near 10s. However, there are significant difference in the magnitude of 
peak SO 3 . The spatial data shows a peak twice as much as the temporal data. The 
wide discrepancy between temporal and spatial results has important implications for 
heterogeneous reactions and for H 2 SO 4 aerosol formulation. This is briefly discussed 
below. 

Note that 503 is produced by the oxidization of SO 2 and then primarily removed 
by heterogeneous reaction on condensed soot particles. As a result, there is an SO% 
mass fraction peak at around 10s. Since the heterogeneous reaction rate is mainly 
related to the thermodynamic condition at the core and almost the same for both 
spatial and temporal simulations, the peak of the S'0 3 mass fraction is determined In- 
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the rate of SO 2 oxidization reaction which is strongly affected by SO 2 mass fraction. As 
Figure 10 and 11 suggest, the spatial data has a higher species concentration around 
the vortex center and thus, a higher production rate of SO 3 than the temporal data. 
These differences between the spatial and temporal simulations imply a necessity for 
spatial simulation in the near field regime. 

The implication of the higher concentration of SO 3 (in the spatial LES) is that 
it could result in an increase in the H 2 SO 4 aerosol formation. Although this feature 
has not been included here, the spatial LES prediction appears to partly explain the 
discrepancy between experimental data and earlier numerical studies. As noted earlier 
[Fahey et ai. 1995], data suggests a 45 times higher H 2 SO 4 aerosols concentration in the 
vortex-plume than the earlier calculations [ Brown et al . . 1996]. The present prediction 
of higher SO 3 concentration in a spatial LES (in contrast to temporal or parabolic space 
marching) seems to suggest that this discrepancy may be partly due to the numerical 
approach. However, to confirm this we need to carry out these simulations, including 
H 2 SO 4 binary nucleation and coagulation aerosol models. This effort is underway and 
will be reported in the near future. 

Some integrated quantities are also analyzed. The mass of entrained NO and SO 3 , 
into the wing-tip vortex are shown in Figure 12. Temporal simulation shows that a 
larger amount of NO and SO 3 is entrained than the spatial simulation does. Since the 
relative reaction rate for NO is small. NO entrainment can be used to represent the 
entrainment process. The difference between the spatial and the temporal simulations 
may be due to the axial velocity deficit of the jet plume. The deviation between the 
two types of simulations begins around x = 12s, where most of outboard engine plume 
is entrained into the wing-tip vortex and the inboard engine plume splits into two 
parts. One part eventually gets entrained into the wing-tip vortex while the second part 
becomes detrained from the wing-tip vortex. The chemical processes that result due to 
the .VO concentration will be different in these two regions. 
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The mass entrainment of SO3 into the wing-tip vortex shows a higher peak for the 
spatial case. The discrepancy in the SO3 concentration in the vortex is particularly 
bothersome as noted earlier. The increased number of aerosol measured in the wake 
[Fahey et «/.. 1995] and the potential increase in the condensed sulfuric acid (due to 
the reaction between SO3 and condensed water) forming on these particles have been 
identified as one process that could impact long term global atmosphere chemical 
balance. The high peak for the spatial case again suggests that simulation boundary 
conditions can affect the predictions especially in the near- field plume. 

In order to study the effect of heterogeneous reactions of HNO3 and H2SO4 on 
plume evolution, temporal and spatial simulations without heterogeneous reactions were 
also carried out. Since the heterogeneous reactions have negligible effect on AO and 
SO2 concentrations, only SO3 data is analyzed. Figure 13 shows the comparison of SO3 
mass fraction at vortex core and entrainment into the vortex core between spatial and 
temporal simulations. Both simulations show much higher SO3 concentration than the 
simulations with heterogeneous reactions. Also, the concentration keeps increasing in 
the downstream direction. The reason is that without heterogeneous reactions ( which 
converts SO3 into H2SO4), SO3 concentration keeps increasing due to the production 
through SO2 oxidization mechanism. Both simulations also show a much higher S' 0 : 3 
entrainment than the simulation with heterogeneous reactions. This is as expected. 
The spatial and temporal results differ much more without heterogeneous kinetics. This 
implies that the spatial simulation yields a higher removal rate of gas phase SO3 inside 
the wing-tip vortex when heterogeneous reactions are included. Furthermore, the result 
suggests that more condensed H 2 SO 4 is formed in the spatial case than the temporal 
simulation. 


Figure 13 



5. Conclusions 

In this study, fluid mixing and chemical reactions in the near-field wake ot B737 
and B747 aircrafts have been investigated. The plume- vortex regime is studied using 
a 3D parallel LES code. Simulations of engine plume interaction with the wing-tip 
vortex of the B737 and B747 have been carried out and comparisons with available 
data have been used to demonstrate the ability of the LES method. Detailed gas phase 
kinetics and simple heterogeneous processes have been included to evaluate the effect 
of turbulent mixing on the chemical processes in the plume-vortex regime. Finally, to 
determine the impact of temporal assumptions (often used in such flow simulations) 
comparisons with spatial simulations were carried out. 

Many features of the computed flow field show reasonable agreement with 
experimental data. The plume-vortex simulation for the B737 shows qualitative 
agreement with LIDAR data and the variation of average SO 2 concentration with 
axial distance for the Bi4/ shows reasonable quantitati\e agieement with axailabh 
data. The entrainment and cletrainment of engine plumes into the wing-tip vortex were 
numericallv captured in the spatial simulation (for the B747) and analyzed. 1 he impact 
of fluid mechanics on the chemical process has been estimated. Results indicate that 
a significant difference between spatial and temporal simulation exists (especially in 
the predicted S0 3 concentration). This difference has an important implication for the 
prediction of sulfuric acid aerosols in the w^ake and may partly explain the discrepancy 
between past numerical studies and the data. Analysis of the results also shows that 
spatial simulation is more suitable for the near field interaction process. For example, 
it was determined that detrainment (observed in the B747 contrail) was only captured 
using the spatial model. 
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Figure Captions 


Figure 1 . Setup of the numerical simulations for the spatial LES simulations. 


Figure 2. Comparison of averaged S0 2 mole fraction decay for different turbulence 
inflow conditions. 


Figure 3. Comparison of B737 scalar concentration with LIDAR measurement (data 
from Poole et al. [1996]). 


Figure 4. Comparison of scalar concentration predicted by spatial and temporal simu- 
lations at a plume age of 5 seconds (30 semispans). 


Figure 5. Comparison of averaged S0 2 mole fraction decay with experimental data. 

Figure 6. Instantaneous vorticity magnitude contours at different times for the spatial 
LES of the B747 wake. The contour interval is 2 sec -1 . 


Figure 7. Vorticity magnitude contours at different locations in the temporal simulation. 
The contour interval is 2 sec -1 . 

Figure 8. Vorticity magnitude contours at different locations in the spatial simulation. 
The contour interval is 2 sec -1 . 



Figure 9. Variation of the detrainment ratio of NO species with the downstream dis- 
tance for B747. 

Figure 10. Maximum NO mass fraction inside and outside the vortex core (vortex core 
has a radius of 4m). 

Figure 11. NO and S'0 3 mass fraction at the vortex core center. 

Figure 12. NO and SO 3 mass entrainment into the vortex core. 

Figure 13. S0 3 mass fraction at the vortex center and entrainment into the vortex core 
with only gas phase chemistry. 



Tables 

Table 1 . Data for B737 ancl B747. 


type 

wing 

maximum 

vortex 

vortex 


span 

circulation 

core radius 

separation 

B737 

30m 

100m 2 /.s 

2m 

26m 

B747 

60m 

600m 2 /s 

4 m 

54m 
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Table 2. Initial conditions for the plume-vortex 
calculations of B747. 

B747 

Plume Ambient condition 


Temperature 

500 1\ 

217K 

velocity! m/sec) 

247 

0.0 

pressure! atm) 

0.17 

0.17 

mole fraction 



CO 

2.4x1 0” 5 

2.0xi0” 8 

C0 2 

3.2x1 0” 2 

3.5xl0" 4 

H 2 

l.OxlO” 6 

5.2x10”” 

Oi 

0.159 

0.2098 

OH 

l.OxlO” 5 

2.8xl0” 13 

h 2 o 

3.0xl0” 2 

4.Sxl0- 6 

Ni 

0.779 

0.789 

NO 

4.3xl0” 5 

5.9xl0” u 

. X0 2 

4.8xl0“ 6 

4.3xlO” 10 

-V 2 0 5 

0.0 

0.0 

HNOz 

0.0 

3.3x10-® 

HONO 

0.0 

0.0 

S0 2 

6.9xl0” 6 

0.0 

so 3 

0.0 

0.0 

0 3 

0.0 

2.5x10”" 
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Figure 1: Setup of the numerical simulations for the spatial LES simulations. 
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Plume age (second) 


Figure 2: Comparison of averaged S0 2 mole fraction decay for different turbulence inflow 
conditions. 


LXDAR M«surtB«nt (B737) Predictions 



Figure 3: Comparison of B737 scalar concentration with LIDAR measurement (data from 
Poole et al. [1996]). 
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Temporal simulation 
Ratio Values 


2.0 +.0 6.0 8.0 10. 


Figure 4: Comparison of scalar concentration predicted by spatial and temporal simula- 
tions at a plume age of 5 seconds (30 semispans). 





Figure 5: Comparison of averaged SO 2 mole fraction decay with experimental data. 



(a) 



Figure 6: Instantaneous vorticity magnitude contours at different times for the spatial 
LES of the B747 wake. The contour interval is 2 sec~ x . 
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(a) x=12s (b) x=24s 


Figure 7: Vorticity magnitude contours at different locations in the temporal simulation. 
The contour interval is 2sec . 



(a) x=12s (b) x=24s 


Figure 8: Vorticity magnitude contours at different locations in the spatial simulation. 
The contour interval is 2 sec -1 . 






Figure 9: Variation of the detrainment ratio of NO species with the downstream distance 
for B747. 


5 .0e-05 



Figure 10: Maximum NO mass fraction inside and outside the vortex core (vortex core 
has a radius of 4m). 
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Figure 11: NO and SO3 mass fraction at the vortex core 
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Figure 12: NO and SO3 mass entrainment into the vortex core. 


S0 3 mass fraction jj S0 3 mass entrainment (kg) 
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Figure 13: S0 3 mass fraction at the vortex center and S0 3 entrainment into the vortex 
core with only gas phase chemistry. 




Numerical Modeling of Chemistry, Turbulent 
Mixing and Aerosol Dynamics in Near-Field 

Aircraft Plumes 

Z. Wang and J.-Y. Chen 
Department of Mechanical Engineering 
University of California 
Berkeley, CA 94720 

A detailed mechanism including gas phase HO x , NO x , and SO x chemistry 
between the aircraft exhaust and the ambient air in near-field aircraft plumes 
is compiled. A reduced mechanism capturing the major chemical pathways 
is developed. Predictions by the reduced mechanism are found to be in 
good agreement with those by the detailed mechanism. With the reduced 
chemistry, the computer CPU time is saved by a factor of more than 3.5 for 
the near-field plume modeling. Distributions of major chemical species are 
obtained and analyzed. The computed sensitivities of major species with 
respect to reaction step are deduced for identification of the dominant gas 
phase kinetic reaction pathways in the jet plume. The potential interaction 
of subgrid turbulent mixing with chemistry is investigated and found to have 
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minor impact on the major species evolution. 

A comprehensive aerosol model is developed and it is coupled with chem- 
ical kinetics and the axisymmetric turbulent jet flow models. The integrated 
model is used to simulate microphysical processes in the near-field jet plume, 
including sulfuric acid and water binary homogeneous nucleation, coagula- 
tion, non-equilibrium heteromolecular condensation, and sulfur-induced soot 
activation. 

The formation and evolution of aerosols are computed and analyzed. The 
computed results show that a large number of ultrafine (0.3-0.6 nm in ra- 
dius) volatile H2SO4-H2O embryos are generated in the near-field plume. 
These embryos further grow in size by self coagulation and condensation. 
Soot particles can be activated by both heterogeneous nucleation and scav- 
enging of H2SO4-H2O aerosols. These activated soot particles can serve as 
water condensation nuclei for contrail formation. Conditions under which 
ice contrails can form behind aircrafts are studied. The sensitivities of the 
threshold temperature for contrail formation with respect to aircraft propul- 
sion efficiency, relative humidity, and ambient pressure are evaluated. The 
computed aerosol properties for different extent of fuel sulfur conversion to 
S(VI) (S 0 3 and H2SO4) in engine are examined and the results are found to 
be sensitive to this conversion fraction. 
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Chapter 1 


Introduction and Literature 
Review 


Research on the impact of aircraft emission on the atmosphere was initiated due to 
interest in commercial supersonic fleet for trans- Atlantic and Pacific rim [Stolarski and 
Wesoky, 1995; WMO, 1995]. Research effort was later extended to subsonic aircraft 
emission impact due to their annual growth rate of 5-7% and environmental concern 
by communities [Schumann, 1995a; Friedl, 1997]. Aircraft emissions in heavily tra- 
versed air corridors may have major atmospheric perturbation. Main concern is the 
effect on O 3 distribution and climate effects. Since there is considerable uncertainty 
about aviation’s impact on the atmosphere, NASA launched the Atmospheric Effects 
of Aviation Project (AEAP), which includes elements of Subsonic Assessment (SASS) 
and Atmospheric Effects of Supersonic Aviation (AESA). The goal of NASA’s SASS 
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project is to develop scientific basis for assessment of atmospheric effects caused by 
engine exhaust emissions, to assist international ozone and climate assessment. The 
assessments serve as necessary recommendations for regulatory agencies, i.e., EPA 
and FAA, as a basis for possible emission standards, and as a guide for technology 
development. 

The emphasis in the assessments of the aviation impact has long been put on the 
emission of nitrogen oxides (NO*) and their potential impact on ozone. Recent in- 
situ emission measurements [Fahey et al., 1995; Schumann et al., 1996; Petzold et al., 
1997] have demonstrated that a large amount of volatile aerosols and a small amount 
of nonvolatile particles are produced in nascent exhaust plumes of subsonic and super- 
sonic aircrafts. Researchers are now especially concerned about the presently largely 
unknown global chemical effects of aircraft-produced particles and soot [Weisensten 
et al., 1996; Bekki, 1997] and how these particles could affect cirrus cloud formation 
[Jenson and Toon, 1997] and thus radiative balance, which in turn could affect the 
climate. 

For simulation purpose, an aircraft wake is divided into three parts: the near-field 
jet plume regime, the vortex interaction regime, and the vortex break-up regime. 
The relationship between different numerical models is displayed schematically in 
Figure 1.1. Plume models are used for tracking chemistry and microphysical processes 
in the near field, while global atmospheric models are intended for estimating the 
potential effects of aviation on the atmosphere. Results from the atmospheric models 



3 


[Bekki and Pyle, 1992; Xue et al. 1994; Weisenstein et al., 1996] indicate potential 
significant atmospheric effects depending, to some extent, on the chemical kinetics and 
microphysical processes in the exhaust plume. The two dimensional global modeling 
study of Weisenstein et al. [1996], for example, predicts a two- fold increase in the 
stratospheric aerosol surface area (under non-volcanic ambient conditions) due to 
sulfuric acid nucleation in supersonic High Speed Civil Transport (HSCT) plumes 
if 10% of the emitted S0 2 is converted to H 2 S0 4 shortly after emission. Near-field 
plume studies will similarly be important for subsonic assessments. 

The present work is focusing on the simulation of near-field jet plume, which 
extends from the engine nozzle exit plane to the location where the plume begins to 
interact with the vortex shed downstream from the wing tip. Near-field modeling 
results serve as an input to the later larger-scale modeling. In the near field, major 
chemical oxidation of gas phase HO x , SO x , NO* by hydroxyl radical OH takes place, 
and volatile sulfate aerosols are generated through binary homogeneous nucleation. 
Chemical activation of soot into condensation nuclei (CN) also occurs in this near-field 
regime. 

A study of near field jet plume with detailed chemical kinetics was previously 
conducted by Miake-Lye et al.[1993]. The flow field was modeled by the JANNAF 
Standard Plume Flowfield model (SPF-II) [Dash et al., 1990] which is an axisymmetric 
parabolic Navier Stokes flow code with a two-equation turbulence model. Karcher 
et al. [1996a] further studied in detail the specific features that may result from 
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assess atmospheric impact 


Figure 1.1: Schematic of relationship among different numerical models in simulating 
aircraft emission impact. 

the coupling between chemistry and turbulent diffusion. More recently, Brown et al. 
[1996a] and Karcher [1996b] investigated the potential importance of microphysical 
processes in the near field. The interaction of subgrid turbulent mixing with chemistry 
was not included in any of the previous simulations. A trajectory model was used 
to study aerosol dynamics by Karcher [1996b] and by Yu and Turco [1996, 1997]. 
Therefore, the flow field is not coupled to the chemistry and microphysical processes 
and no radial resolution is offered. In addition, Yu and Turco’s [1996, 1997] work did 
not include chemical kinetics. 

In the wake- vortex regime, the flow is highly three-dimensional. Large-eddy simu- 
lation (LES) is believed necessary to accurately resolve the complicated fluid dynam- 
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ics. However, LES is computationally expensive and it is not feasible to incorporate 
detailed chemistry. A reduced chemistry is therefore required to alleviate the compu- 
tational burden. 

The objectives of the present work are: 1) to compile a detailed chemical mecha- 
nism for modeling reactions between engine exhaust and ambient air; 2) to develop a 
reduced chemistry for LES in the wake-vortex regime; 3) to study chemical transfor- 
mation in the near field and to predict distribution of key species; 4) to explore the 
interaction between turbulence and chemistry using the probability density function 
(pdf) approach; 5) to develop a comprehensive aerosol model for studying microphys- 
ical processes of aerosols in the near field plume; 6) to investigate the sensitivity of 
particle properties to less-known parameters; and 7) to study the threshold condition 
for contrail formation behind aircrafts and to access its sensitivity. The present mod- 
eling provides detailed predictions of gas phase species as well as aerosol particles, 
which serve as an input for the follow-on model. 
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Chapter 2 

Model Formulation 

2.1 Model Structure 

Figure 2.1 shows schematically the chemical and microphysical transformation 
processes of aerosols modeled in the present work. SOx and soot particles are emitted 
from the nozzle exits of jet engines at high temperatures (500 K - 600 K). Sulfuric acid 
H 2 SO 4 is generated through gas phase chemical reactions of SOx with OH and H 2 0. 
Soot particles can be chemically activated by adsorption of SO 3 and H 2 SO 4 , leading 
to the formation of a partial or full liquid H 2 S0 4 -H 2 0 coating. As the plume cools 
by the entrainment of the cold ambient air, volatile liquid H 2 SO 4 -H 2 O aerosols are 
formed through binary homogeneous nucleation. These volatile aerosols grow in size 
through either condensation of gaseous H 2 SO 4 and H 2 0 or by self coagulation. The 
volatile aerosols can interact with soot particles through coagulation, which provides 
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Figure 2 . 1 : Schematic of particle microphysical processes in near-field aircraft plumes 
as simulated by the aerosol model. 


another pathway for the sulfur induced soot activation. H2SO4 and H2O can condense 
onto the coated part of soot particles. 


2.2 Model Input 

To initialize our numerical calculations and to convert emission indices of primary 
exhaust products into mole fractions, we need engine exit-plane properties of aircrafts. 
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Table 2.1: Exit plane parameters for B747 engine 



Exhaust Core Flow 

Bypass Air Flow 

Ambient 

Fuel flow rate mp, kg/s 

0.795 



Exit radius, m 

0.44 

0.435 


Velocity, m/s 

475.7 

316.3 

237 

Temperature, K 

547.3 

253.4 

219.2 

Pressure, atm 

0.2361 

0.2361 

0.2361 


We will concentrate on the CF6-80 C2B1F engine used by the Boeing 747 aircrafts. 
The formation and evolution of chemical species and particles are expected to develop 
similarly behind other engines, with slight changes in the length scales due to different 
nozzle geometry and axial exhaust velocity. The results of this work are considered 
to be representative of other types of aircrafts qualitatively. 

Recently, Schumann [1995b] has estimated the relevant nozzle exit plane parame- 
ters for the B747 engine in such a way that they are consistent with the integral mass, 
momentum, and energy balances of the jet core and bypass flows. Table 2.1 lists the 
parameters adapted from Schumann [1995b] for a B747 at the cruising condition of 
P a = 0.2361 atm and T a = 219.2 K with airspeed u a = 237 m/s. 

Given the emission indices Elk of the exhaust product k in gram per kilogram of 
fuel burnt, the species volume mixing ratios (mole fractions) at the nozzle exit plane 
can be derived from the following relation based on mass conservation 


%k, 0 


Eh rhp 1 
m k V N a ' 


( 2 . 1 ) 


where m* is the molecular mass of species k, mp is the engine fuel flow rate, V = 
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A COT efl u corcfi is the volumetric flow rate of the core flow at the exit plane. The subscript 
0 denotes the engine exit plane, and N a = P a f{k B Tcarefi) is the air number density at 
engine exit plane. With the data in Table 2.1, we find for a B747, 


x fc ,o = 5.27 x 10' 4 


Eh [g/kg] 
M fc [g/mol] ’ 


( 2 . 2 ) 


where M* is the molar mass of species k. 

In Table 2.2, we present the initial mole fractions at the exhaust core based on 

equation (2.2). The mole fractions in the ambient air are prescribed as the boundary 

conditions in the numerical model. These ambient compositions are extracted from a 

mesoscale chemical transport model based on the RADM mechanism [Chang et al., 

1987; Petry et al., 1994] and represent a region with dense air traffic, e.g., at northern 

midlatitude (50°N) summer conditions at noon. The choice of the boundary values 

is not crucial since we only focus on the short-term evolution of the near-field jet 

regime. The ambient water mole fraction listed in Table 2.2 corresponds to a relative 

humidity of RH = x J^tO,aPa ^ w h ere the subscript a denotes the ambient 
E fat, Hi O \Ea ) 

condition and the superscript oo denotes a flat surface. 

The emission indices of water vapor and carbon dioxide can be estimated relatively 
accurately assuming stoichiometric jet fuel combustion. Jet fuel sulfur contents (FSC) 
typically vary from 0.1 to 1 g/kg [Busen and Schumann, 1995], with a broad average 
around 0.4g/kg (400ppmm) and a maximum specification limit of 3 g/kg [American 
Society for Testing and Materials (ASTM), 1994]. Assuming that fuel sulfur is totally 
oxidized to sulfur dioxide during combustion in the engine, the SO 2 emission index 



Table 2.2: Initial and boundary conditions of gas species for B747 



Exhaust Core Flow 

Ambient Air 

Mole Fraction 
0 

0 

1.1 x 10" 14 

o 2 

0.135 

0.21 

0 3 

0 

3.0 x 10- 7 

H 

0 

0 

h 2 

0 

9.0 x 10“ 7 

OH 

6.2 x 10~ 6 

2.8 x 10" 13 

ho 2 

0 

3.6 x 10- 12 

h 2 o 

3.6 x 10" 2 

8.8 x 10~ 5 

h 2 o 2 

0 

2.4 x 10" 11 

NO 

1.3 x 10“ 4 

8.1 x 10" 11 

no 2 

6.9 x 10" 6 

6.5 x 10“ u 

N0 3 

0 

8.9 x 10“ 15 

n 2 o 5 

0 

1.0 x 10~ n 

hno 2 

0 

3.4 x 10" 13 

hno 3 

0 

2.0 x 10" 9 

CO 

9.4 x 10“ 5 

4.0 x 10- 8 

co 2 

3.8 x 10" 2 

3.3 x 10" 4 

so 

0 

0 

so 2 

5.8 x 10~ 6 

9.1 x 10- 12 

so 3 

0 

0 

hso 3 

0 

0 

h 2 so 4 

0 

0 



11 


can be directly calculated with EIso 2 = FSC • Mso 2 /^s- In this work, we use 
EIso 2 = 0.7 gj kg. A sensitivity study with different levels of S0 3 emission will be 
presented in Chapter 4. The emission indices of the other two major air pollutants 
resulting from combustion, CO and NO* (= NO + NO 2 ), are also less well known 
because they depend on several factors including combustion temperature, fuel-air 
ratio, hydrocarbon oxidation process, and fuel nitrogen content. In the present work, 
we employ EIco — §9 /kg and EI^ox equivalent to 12 g of NO x per kg fuel burnt. 
NO z emission is assumed to consist of 95% NO and 5% NO 2 on a molar basis. 

The hydroxyl radical OH is a key gaseous species driving major radical chemical 
reactions in the cooling jet plume. Unfortunately its emission index is currently only 
known with considerable uncertainties from simulations of the internal engine flow 
and combustion [Miake-Lye et al., 1993; Stolarski and Wesoky, 1995]. Karcher [1996] 
showed that EI 0 h = 0.2 g/kg fuel, corresponding to 6.2 ppmv from equation (2.2), 
can be used as a reasonable reference case. The abundances of all other radicals at 
the exit plane are set to zero mainly because their emission indices are likely to be 
insignificant or simply unknown. 

Typical size distributions of nascent soot particles are assumed to be log-normal, 
with mean radii in the range of 15-30 nm, modal widths of 1. 4-1.6, and number 
densities of 10 6 -10 7 cm -3 [Hagen et al., 1992; Petzold and Schroder, 1997]. In this 
work, we take the median radius as 20nm, the modal width as 1.4, and the total 
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number density as 1.6 x 10 6 cm* 3 . Since the number density can be expressed as 

\t EIaoot™-F (n o\ 

Naoot = *77 > V 2 - 3 ) 

ym 30 ot 

where the mass of a soot can be approximated as m aoot = p SO ot * 4/37rf 3 with p soot ~ 2 
g/cm 3 for graphitic carbon. The soot emission index is estimated to be 0.04 g of soot 
per kg of fuel. This corresponds to an emission index of El soo t/{p SO ot * 4/3tt f 3 ) = 
0.6 x 10 15 soot particles per kg of fuel. 

2.3 Solution Method 

To simulate the near-field jet plume, we assume the flow field is axisymmetric, 
steady-state, and parabolic. Isobaric expansion at a prescribed ambient pressure P a 
is assumed. The turbulent mean velocity field is modeled by the Reynolds average 
Navier Stokes equations with a second moment turbulence closure model, which di- 
rectly solves the modeled transport equations for the Reynolds stresses. A marching 
downstream algorithm is used and the details can be found in Chen et al. [1987, 
1989]. When the turbulence interaction with chemistry in a subgrid is neglected, the 
mean concentrations and temperature can be solved by the conservation equations 
of mean species and energy. The mean chemical reaction rate can be expressed in 
terms of mean temperature and mean concentrations u> k = f(f,Ci). To account for 
the subgrid interaction of turbulence with chemistry, an additional term of subgrid 
molecular mixing is needed, and the mean chemical reaction rate cannot be expressed 
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in general as a function of mean temperature and mean species concentrations. An 
efficient way to resolve this problem is by using the probability density function (pdf) 
approach. In general, the instantaneous composition of a reacting mixture of gases 
can be determined from the set of kk + 1 scalars comprising mass fractions and 
enthalpy fax,t) = fa, fa, ••• ,fak,fak+u where kk is the total number of gaseous 
species. With %p = ,ipkk+ 1 being the composition space corresponding to 

<t>, the mean value of any quantity (Q(fa) can be determined from the probability 
density function of joint scalar fa Pj,(ip;x,t), by (Q(fa) = f Pj(iP)Q (fa dip, where 
P^{fax,t) is defined as P$(fax,t)dip = Probability^ < 4> < ip + dPp). In variable- 
density turbulent flows, the density-weighted average of a function Q is defined as 
Q = ( pQ)/(p ) = f P+(i>)p{ip)Q{fadip/(p) = fPi$)Q$)dfa Therefore, the density- 
averaged reaction rate can be expressed as u k = f Pj(ip)uk(ip)dip. 

The density-weighted joint scalar pdf can be solved by the following transport 
equation [Pope, 1985]: 


(P)uk 


dP^jfax) 

dx k 


^ ft r _ _ T 

- E ^ 

a=l 


_ _d_ 

dx k 

** + l 

- E 


[WK|« = ^)P,W;f)' 


(2.4) 


a 2 


[0>)(£«<>l^ =$)Pi(t \ *)] . 


where () denotes the mean value, e Q0 is the scalar dissipation rate given by e a0 = 


r dfadfa wherg r . g the molecu i ar diffusivity, (e Q0 \<p = ip) represents the average 
dx k dx k 
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scalar dissipation rate conditioned on a certain set of scalars, and S Q is the chemical 
source term for scalar ip a . The first two terms on the right hand side represent the 
effects of chemical reaction and turbulent transport on the joint scalar pdf, respec- 
tively. Both of these terms can be calculated directly with no need for modeling. The 
last term stands for the effect of microscale mixing on the pdf which requires closure 


approximations. In this study, the effect of microscale mixing on the pdf is modeled 
by the modified Curl’s mixing model [Janika et al., 1979], 


- V a A " \(eafi\* = $)Pi(fr2 


^mix t Jib* J \p n J 


where H is the transitional probability defined as 


H(xP', <|V») = rrr, ~~ ~77 T for € [V>', i>"] 

\W -V\ (2.6) 

= 0 else wise. 


The amount of mixing is modeled through a characteristic mixing time which is 

related to the turbulence time scale for micromixing, r = k/i, where k = 0.5 (t^u*) 

du" du" 

is the turbulence kinetic energy and i = ^(-^- L -x- L ) is the turbulence dissipation 

OX/c OXk 

rate with v being the kinematic viscosity. Reactive flows of practical interest usually 
involve many species, leading to a large dimensionality of P(ip\,rp 2 , • • ■ , i>kk+i)- The 
computational expense of finite difference method is found to rise exponentially with 
the dimensionality and it is impractical to solve a variable with dimensionality greater 
than three. However, the computational expense of the Monte Carlo technique rises 
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only linearly with the dimensionality, which is the best that can be achieved by 
any algorithm. Therefore, the Monte Carlo technique is employed to solve the pdf 
evolution equation (2.4). 

The turbulent transport term can be modeled by a gradient diffusion model using 
turbulence time and fluctuating velocities from the Reynolds stress model, 

<t*w = = - c 4^ SP dx e ' x) ’ (2 - 7) 

where C 3 is an adjustable constant (= 0.25 in jet flows), u'j. = Uk — Uk is the density- 
weighted velocity fluctuation and (u" k \<j> = xj>) is the average fluctuation conditioned 
on a certain set of scalars. 
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Chapter 3 

Gas Phase Chemistry 

3.1 Detailed Chemistry 

3.1.1 Compilation of Detailed Chemical Reaction Mechanism 

Table 3.1 and 3.2 list the bimolecular and termolecular gas phase reactions em- 
ployed in our model, which include 60 chemical reactions and 24 species. The reaction 
kinetic parameters are selected to cover the temperature range from 1000 K to 200 
K as required by the jet plume temperature history. Photochemical reactions are 
excluded as the time scale of the fastest photolytic reaction, NO 2 + hv — t NO + O, 
is about 2 minutes which is larger than the plume age of just a few seconds. Hy- 
drocarbon reactions are also excluded because the oxidation potential of the young 
exhaust plume is too small to affect the emitted hydrocarbons. We also do not at- 
tempt to model ion chemistry that may influence the radical balance in the very early 
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jet regime. 

Most of the reaction rate parameters are extracted from the chemical kinetics 
database collected by the National Institute of Standards and Technology (NIST) 
[Mallard et al., 1994]. Rate constants for the reaction of HNO 2 with H and O (reac- 
tions (R28) and (R29)) are taken from data by Tsang et al. [1991]. JPL database 
[Demore et al., 1992] on chemical reaction rate constants provides accurate rate data 
only for common atmospheric temperature below ~ 300 K. To model the kinetics of 
the reactive jet flow properly, however, we need to employ rate constants which are 
also valid at higher temperatures. In some cases, a simple extrapolation of the JPL 
rates would give similar answers, but this may not hold for all reactions. 

3.1.2 Chemical Kinetics Package 

The Chemkin package [Kee et al., 1989] is a software tool designed to facilitate 
simulations of elementary chemical kinetics. An important advantage of the general- 
purpose and problem-independent structure of Chemkin is that it allows us to work 
with the same chemical input for different application codes. As indicated in Fig- 
ure 3.1, this package consists of two major FORTRAN codes, the Interpreter and 
the Gas-Phase Subroutine Library, along with the Thermodynamic Database. When 
the Interpreter is executed, it reads in the user’s symbolic description of the reac- 
tion mechanism and then extracts the appropriate thermodynamic information for 
the species involved from the Thermodynamic Database. The major output of the 
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Table 3.1: Gas Phase Chemical Reaction Mechanism: Binary Reactions 


Reaction 


Rate Constant (cm 3 molecule 1 sec 1 ) 


R1 

R2 


O + O3 
H + 0 3 


20 2 

OH + O 2 


1.21 x 10 -1I exp(-2125/T) 

1.15 x 10~ 10 exp(-436/T) 

8.10 x 10~ 21 T 2 8 exp( — 1 950/T) 
2.8 x 10 -1 °exp(-440/T) 


R3 H + OH — 4 O + H 2 
R4 H + H0 2 — 4 OH + OH 
R5 H + H0 2 — 4 H 2 + 0 2 
R6 H + H0 2 —4 H 2 0 + O 
R7 OH + O — 4 H + 0 2 
R8 OH + 0 3 — 4 H0 2 + 0 2 
R9 OH + H 2 — 4 H 2 0 + H 
RIO OH + OH —4 H 2 0 + O 
Rll OH + H0 2 — 4 H 2 0 + 0 2 
R 12 OH + H 2 0 2 —4 H 2 0 + H0 2 
R13 H0 2 + O — 4 OH + 0 2 
R14 H0 2 + 0 3 — 4 OH + 20 2 
R15 HO 2 + HO 2 — 4 H 2 O 2 + O 2 
R16 H 2 0 2 + O — 4 OH + H0 2 
R17 H 2 0 2 + H — 4 OH + H 2 0 
R18 H 2 0 2 + H — 4 H0 2 + H 2 
R19 NO + 0 3 — 4 N0 2 + 0 2 
R20 NO + H0 2 —4 N0 2 + OH 
R 21 NO + N0 3 —4 N0 2 + N0 2 
R22 N0 2 + O — 4 NO + 0 2 
R23 N0 2 + 0 3 —4 N0 3 + 0 2 
R24 H + N0 2 —4 OH + NO 
R25 N0 2 + N0 3 — 4 NO + N0 2 + 0 2 
R26 N0 3 + 0—4 N0 2 + 0 2 
R27 OH + N0 3 —4 H0 2 + N0 2 
R28 HN0 2 + O — 4 OH + N0 2 
R29 HN0 2 + H — 4 N0 2 + H 2 
R30 HN0 2 + OH — 4 H 2 0 + N0 2 
R31 HN0 3 + O — 4 OH + N0 3 
R32 HN0 3 + OH — 4 H 2 0 + N0 3 
R33 SO + 0 2 — 4 S0 2 + 0 
R34 SO + O 3 — 4 SO 2 + O 2 
R35 SO + OH —4 S0 2 + H 
R36 SO + N0 2 —4 S0 2 + NO 
R37 S0 2 + 0 3 — 4 0 2 + S0 3 
R38 S0 3 + O — 4 0 2 + S0 2 
R39 S0 3 + H 2 0 —4 H 2 SO 4 

R40 HS0 3 + 0 2 — 4 H0 2 + S0 3 
R41 CO + OH —4 C0 2 + H 


6.9 x 10 -n exp(-636.9/T) 

2.8 x 10~ 12 T 046 exp(— 677.9/T) 
1.83 x 10 -11 exp(173.3/T) 

1.9 x 10 _12 exp(-1000/T) 

1.11 x 10 _16 T 164 exp(— 1589/T) 
8.34 x 10- 17 T 1M exp(355/T) 

5.09 x 10 _11 exp(72.6/T) 

2.13 x 10 -13 T 0-47 exp(— 179.8/T) 
2.71 x 10 _11 exp(-224/T) 

1.4 x 10 _14 exp(-600/T) 

2.2 x 10 _13 exp(600/T) 

2.33 x 10 _11 exp(-2814/T) 

1.7 x 10~ n exp(-1800/T) 

1.77 x 10“ 1 1 exp(-2890 /T) 

2.14 x 10 _12 exp(-1408/T) 

3.7 x 10~ 12 exp(240/T) 

1.8 x 10 _u exp(110/T) 

6.5 x 10 -12 exp(120/T) 

1.2 x 10 _13 exp(-2450/T) 

1.4 x 10“ 10 

1.91 x 10 _13 exp(-1696/T) 

1.0 x 10~ n 

2.3 x 10- 11 

2.0 x 10 -11 exp(-3000/T) 

2.0 x 10“ n exp(-3700/T) 

1.8 x 10" 11 exp(-390/T) 

3.0 x 10" 17 

4.02 x 10 -14 exp(317.7/T) 

1.55 x 10 _13 exp(-2288/T) 

4.3 x 10 -12 exp(-1148/T) 

8.59 x 10“ n 

1.4 x 10" 11 

3.0 x 10~ 12 exp(-7000/T) 

3.17 x 10“ 1 1 exp(-4455 /T) 

1.2 x 10" 15 

1.23 x 10 _12 exp(-316.8/T) 

4.44 x lQ- 16 T°- 98 exp(94/T) 
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Table 3.2: Gas Phase Chemical Reaction mechanism: Three-Body Combination and 
Thermal Decombination Reactions 


Reaction Low-Pressure Limit k 0 High-Pressure Limit hoc 

R42 O + O Oj 5.2 x 10- 3s exp(900/T) 

5.51 x 10-*'T- 1 3 exp(152.6/T) 
i.3o x io- 29 r- ! 


R43 O + 0 2 0 3 

R44 H + O -^+ OH 

R45 H + Oj A H0 2 
R46 H + H H a 

R47 H + 0H-^+H 2 0 
R48 OH + OH H 2 0 2 + 0 2 
R49 H0 2 + H0 2 H 2 0 2 

R50 NO + O N0 2 

R51 OH + NO A HN0 2 
R52 N0 2 + O NOa 

R53 OH + N0 2 HN0 3 
R54 N0 2 + N0 3 N 2 O s 

R55 N 2 O s N0 2 + N0 3 

R56 HN0 2 OH + NO 

R57 HN0 3 OH + N0 2 

R58 S0 2 + O S0 3 

R59 S0 2 + OH -^+ HS0 3 

R60 CO + O A CP 2 


1.05 x 10” 30 T” 0 7 exp( 144 .2/T) 
1.28x 10 -29 T -1 13 exp(-100.6/T) 
1.538x10- 28 T-> 2I exp(295.3/T) 
6.58 x io- 29 T"° 8 

1.9 x 10" 33 exp(980/T) 

9.09 x 10- 28 T-’ 6 

1.91 x 10- 26 T-’ 9 exp(164.8/T) 
6.21 x 10 -23 T -3 - 22 exp(-570.5/T) 
3.50 x 10- 16 T- S J exp(-840/T) 
6.98 x io- 22 T- 3 4 

1.5 x 10 _6 exp(-9104/T) 

5.0 x 10 6 T” 3 8 exp(-25340/T) 

1.2 x 10 _6 exp(-22270/T) 

1.49 x 10 -32 exp(-601/T) 

1.97 x 10 _32 exp(867.3/T) 

2.49 x 10~ 33 exp(-1550/T) 


4.2 x 10 -12 exp(-183/T) 


7.5 x 10" 11 

2.9 x 10 -lo exp(5821/T) 

2.69 x 10 -1 °exp(-75/T) 

1.5 x 10" 11 

2.2 x 10~ 13 exp(600/T) 

5.43 x 10- 12 T° 3 

5.1 x 10 -12 exp(159.8/T) 

2.65 x 10 _I1 exp(-25.78/T) 

1.84 x 10~ 12 exp(587.2/T) 

6.40 x 10-' 3 T° 2 

5.9 x 10 n exp(-9370/T) 

1.20x 10 19 T" 1 23 exp(-25010/T) 
2.14 x 10 15 exp(-24270/T) 

1.19 x 10~ 9 
2.0 x ltr 12 

2.66 x 10~ 14 exp(-1459/T) 


The effective second-order reaction rate constant is [Demore et al., 1992]: 


■ t) =(t 


fco 




+ [log 10 (*0C7’}[M]/fc oo {T))] 2 >- 


For third-order low and Becond-order high pressure limits, the units are cm 6 mo]ecule“ 2 sec~ 1 and 
cm 3 molecuIe‘ - ^ec" 1 , respectively. For second-order low and first-order high pressure limits, the units are cm 3 sec“ 1 
and sec" 1 , respectively. 


Interpreter is a linking file, which contains all the pertinent information on the ele- 
ments, species, and reactions of the specific mechanism. The Gas-Phase Subroutine 
Library is a collection of more than 100 subroutines that return information on ele- 
ments, species, reactions, equations of state, thermodynamic properties, and chemical 


reaction rates. 
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Chemkin Package 



Figure 3.1: Flow chart of using the CARM code with Chemkin and PSR code to 
develop reduced mechanisms. 
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3.2 Reduced Chemistry 

In the intermediate regime of the wake, the jet plume will be entrained into the 
wing tip vortex. To model the impact of the wing tip vortex on chemistry in this 
regime, large eddy simulation (LES) may be required in order to model the three- 
dimensional vortex properly. Solving the flow field by LES would demand significant 
computer time, making it impossible to incorporate a detailed chemical mechanism. 
A reduced mechanism capable of describing the key reactions is necessary. 

3.2.1 Perfectly-Stirred Reactor (PSR) Model 

Perfectly-Stirred Reactor (PSR) model will be used in the current study to develop 
reduced chemistry. PSR is an idealized reactor which is assumed to be homogeneous 
(spatially uniform) and well-mixed at the molecular scale. As shown in Figure 3.2, 
there is a continuous mass flow in and out of the reactor as shown in Figure 3.2. 
When the incoming fluids enter the reactor, it is assumed that the fluids are perfectly 
mixed with the fluids existed in the reactor immediately. This assumption is justified 
when the characteristic time of turbulent mixing at both macroscale and microscale 
is much smaller than that of chemistry. Consequently, the rate of chemical process is 
controlled solely by chemical kinetics. Turbulent effects are not accounted for in this 
model because of the assumption of infinitely fast mixing. 

The conservation equations for reactive species at the steady state PSR is given 
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by 


m{Yk-Y;)-u k M k V = 0, (3.1) 

and the energy conservation equation is described by 

KK 

m'52(Y k h k -Y;h' k ) + Q = 0. (3.2) 

*=i 

In these two equations, Y k is the mass fraction of the fcth species (there are a total 
of KK species), M k is the molecular weight of the A;th species, V is the reactor 
volume, m is the mass flow rate through the reactor, u k is the molar production rate 
of A:th species per unit volume by chemical reactions, h k is the specific enthalpy of 
the kth species, and Q is the reactor heat loss. The superscript asterisk indicates the 
inlet conditions. The characteristic flow time is residence time r = pV/m. A hybrid 
Newton/time-integration algorithm is used to solve for the solution and details are 
described by Glarborg et al. [1986]. 


3.2.2 Derivation of reduced chemistry 


A reduced chemistry is feasible if certain intermediate species reach the quasi- 
steady state when their production rates are nearly equal to their consumption rates. 
To identify a species in the quasi-steady state, we apply the criterion: 






max(|^|, |wj|) 


< 1 , 


(3.3) 


where u> k and u k are the production and consumption rate of species it, respectively, 
and X k is the mole fraction of species k. It is noteworthy that due to the rapid 
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temperature changes in the near field of jet plume, a species that does not satisfy 
the quasi-steady state criterion in the near field may reach it in the far field, such 
as radicals OH and SO 3 in the present case. One may construct different reduced 
mechanisms suitable for different wake regimes. 

The present development of reduced chemistry has been carried out using the PSR 
with an interactive Computer-Assisted Reduced Mechanism (CARM) code [Chang, 
1995] for the near field jet plume. Figure 3.1 shows the flow chart of using CARM 
with Chemkin and PSR codes to develop a reduced mechanism. 

Assuming the quasi-steady state for 0, H, H 2 , H 2 O 2 , H0 2 , SO, HSO 3 , NO 3 , and 
N 2 O 5 (9 species out of a total of 24 species), we obtain the following eight-step reduced 
mechanism: 

(G-Rl) 20 3 = 30 2 

(G-R 2 ) 0 2 + 20H = H 2 0 + 0 3 

(G-R3) NO + 0 3 = 0 2 + N0 2 

(G-R3) 20 2 + NO + HN0 3 = OH + 2N0 2 + 0 3 

(G-R 6 ) 20 2 + HN0 2 = OH + N0 2 + 0 3 

(G-R5) O 3 + SO 2 = 0 2 + SO 3 

(G-R7) H 2 0 + S0 3 = H 2 S0 4 

(G-R 8 ) O 3 ■+■ CO = 0 2 + C0 2 

The rates for the global reactions are expressed in terms of elementary reaction rates 
as: 
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UC-R\ — w 3 — + <^e + U>10 — U>13 — 0^16 — b>n — ^20 + w 25 + w 27 + w 29 + ^30 + w 32 

— U) 42 ~ W 43 — U 44 + i0 4S - U> 5 0 ~ 0 > 5 i — 0>52 ~ ^53 + ^56 + <*>57 ~ w 58 ~ ^60 
U>C— R2 = ^6 4" ^9 4“ <*>10 — tc^l] 4" ^12 4" ^17 4" ^30 ^32 4" U} 47 

Uq-RZ = Wi94-a>20 — ^224 _ ^23 — ^24 — 2^25 — ^26 — ^27 _ ‘^364'^504"^5l4-a;524-u;53 — CJ56 — CJ57 

^G-fi4 = W 21 “ ' t^23 4" <^25 4" ^26 4" ^27 — ^52 — ^53 4" w 57 

U>G-R5 — ^28 + ^>29 + ^30 — W 51 4- ^56 

U>G-R6 = ^37 — ^38 4" w 40 4" ^58 

WG-R7 — w 39 

UlQ-RS = W41 4- W60 

Although the species in the quasi-steady state do not appear explicitly in the 
global reactions, their species concentrations are needed in computing the elementary 
reaction rates which contribute to the global reaction rates as shown above. The 
concentrations of quasi-steady state species are computed using the following expres- 
sions: 


[HS0 3 ] = 


MggHggg] 

£ 40 (^ 2 ] 


k 3 i[0][HN0 3 ] ± MjVggs] + fc52[0][^Q 2 ] 4- fcatJVQaHQs] 4- k 32 [0H][HN0 3 ] 
[iV03] “ A: 26 [0] + 4- k„[OH] 4- k b4 [N0 2 ] 4- k 2 1 [NO] 

jb 3 5 [OH][SO] 4- k 7 [OH][Q] 4- fc 9 [OH][i/ 2 ] 4- k A1 [OH][CO] 

[H] - J 

A = k 46 [H ] 4- k 3 [OH] + k 44 {0][M] 4- k lB [H 2 0 2 ] + k 29 [HN0 2 ] 4- k l7 [H 2 0 2 ] 4- k 6 [H0 2 ) 


+ k 5 [H0 2 ] 4- k 47 \OH] 4- k 4 [H0 2 ] + k 2 \0 3 ] + k 24 [N0 2 ] + k 45 [0 2 ) 
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[01 _ MOgUgO] + h\H][OH] + * 6 (ff](ff0 2 ] + k l0 [OH][OH] 

1 J B 

B = * 6 ,33 [S0 2 ] + *38 [S0 3 ] + *44 [H] [M] + * 42 [0] [M] + * 16 [H 2 0 2 ] + *60 [CO] + * 31 [ff NOa] 
+ *2s[-^-^ ^ 2 ] + *26[^V ^ 3 ] + + *58 + * 1 ( 03 ] + k-j\OH\ + *52(-N0 2 ] 

+ k b0 [NO] + *22(-N0 2 ] + *43 [^ 2 ] 

lso1 t»3|Q|[SQ2) 

*34(^3] + * 35 [ 0 //] + *36(^02] + *33(02] 

ri7 , _ k 3l [0]\HN0 3 ] + ku{0][HN0 3 ] + k 31 [0][HN0 3 ] + * 31 [0][tf JV0 3 ] 

^ k 9 [OH] 

lHO]= k^[HQ2\[H02\ ± fcisjgOgHg^] + ku[OH][OH] 

*18 [-W] + *i6[0] + kn[H] + k\ 2 [OH] 

I ho 2 ) = | 

A = k m \H}\H 2 O 2 \+k li \0}\H 2 0 2 \+k v \0H}\N0 1 } + k n [0H]\H. t 0 2 ] + k t \OH}[O i \ 
+ k 40 [O 2 ][HSO 3 ] + k^HWO^ 

B = kf,[H] + k$[H] + *4 [-*f] + * 13 ( 0 ] 4- k 49 [H0 2 ] + ki$[H0 2 ] + * 14 ( 03 ] + k n [OH] 
+ k 20 [NO] 


|W2 o sl = 

*55 

Consequently, the concentrations of the quasi-steady state species can be expressed 
in terms of major species by a set of algebraic equations. This saves computing time 
as there is no need to solve ordinary differential equations (ODE) for these species. In 
addition, the stiffness of the ODEs with the global chemistry derived is reduced since 
the species in the quasi-steady state are often fast-reacting. The Newton iteration 
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m Residence Time, x 

Y*(k) ► Temperature, T 

h*(k) Enthalpy, h(k) 



m 

Y(k) 

h(k) 


Figure 3.2: Schematic of a Perfectly-Stirred Reactor (PSR). The characteristics of the 
reactor is shown, where rn is the mass flow rate, k denotes species k, and superscript 
* indicates the inlet condition. The key parameter is the residence time, r = pV/m. 


is used to solve the coupled nonlinear algebraic equations for concentrations of the 
quasi-steady state species. 


3.3 Sensitivity Analysis 


Sensitivity information is useful for identifying the controlling chemical steps that 
have strong influence on the creation/destruction of a certain species. The sensitivity 
coefficients are displayed as 



(3.4) 


where is the sensitivity of change in the mole fraction of the kth species, X k , due 
to a small change in temperature-independent pre-factor, A { , of the i-th reaction rate 
ki = AiT n exp (~E a /T). 
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The transient well-mixed reactor (WMR) model enables one to study the time 
evolution of chemical kinetics of the engine exhaust gas subject to injection of ambient 
air. A WMR model is employed here to investigate the sensitivity of species with 
respect to reaction step. A schematic diagram for the WMR is sketched in Figure 3.3. 
Similar to the PSR, a global residence time can be defined as 


m R 

tr = 

m 0 


(3.5) 


where m R is the mass of gas inside the reactor and m a is the mass entrainment rate. 
Using tr as a parameter, conservation equations of mass, species, and energy can be 


written as 


drriR ttir 

tr 


dt 

dYni = T(y OJ - y S j ) + 


dt 


dTi 


tr 

N 


URj 

Pr 


Mj, j = l 


(3.6) 

(3.7) 


dt 


" = [(a. . - h n ,)n,] - y, (3 . 8) 

> CpRTRfr', J 3 ° t-T PRCPR 


j = l 


CpRTTlR 


where the subscript R represents quantities inside the reactor and those entrained are 
denoted by the subscript a. If the temperature is fixed, the WMR model is similar 
to the Lagrangian box model [Miake-Lye et al., 1993b; Karol et al., 1994; Karcher, 
1995] when the following relation is used to bridge the two models 


m 0 


m R 


m 


(3.9) 


A is the cross-section area of a plume and K(t) is the growth rate of the plume 
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Figure 3.3: Schematic diagram of a transient well-stirred reactor (WMR) showing 
ambient air entrained into the reactor with a rate of m a . ttir is the mass of gas in 
the reactor, and it increases with time. 


cross-section area in Lagrangian box model equation 

^ = P, - L, - ^(Ci - C“) + fl,, (3.10) 

where C* and C“ are the ith species concentrations inside the box and in the air, 
respectively, Pi and Li are the production and destruction rate by chemical reactions, 
and Ri is the source or sink due to other mechanisms. 
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3.4 Results and Discussions 

The WMR model is used to provide sensitivity information of a certain species 
with respect to various reactions. The entrainment rate of ambient air is prescribed 
from centerline mass increment rate predicted by the pdf jet plume model. Therefore, 
the sensitivity information is informative for understanding qualitatively the center- 
line species distribution in a jet plume. Figure 3.4 shows the computed sensitivities 
of hydroxyl radical, OH, with respect to various reaction steps. The most important 
chemistry for the OH radical is its reaction with NO via the three-body step (R51) 
OH + NO + M — ► HNO2 + M. The generated HNO2 then further reacts with OH 
to produce NO2 by step (R30) HNO2 + OH -4 H 2 0 -I- NO2. The self combination re- 
action step (RIO) OH -I- OH — > H 2 0 -I- O is the next important pathway contributing 
to the OH consumption in the early jet regime when the OH level is high. (R41) 
CO + OH -4 CO2 4- H and (R48) OH 4- OH 4- M -4 H2O2 4- M also play important 
roles in this early stage. The magnitudes of the OH sensitivity coefficient start to 
decline when the OH abundance begins to decrease dramatically. They become very 
small when OH decreases to a very low level. Compared to the major OH consump- 
tion pathways mentioned above, a very small amount of OH is consumed by SO2 via 
step (R59) SO2 4- OH 4- M — ► H2SO3 4- M, which has a very small magnitude of 
OH sensitivity. This explains why the conversion rate of SO2 to H2SO4 is less than 
1-2 % as will be shown later in this section. A small amount of OH can be gener- 
ated through (R20) NO -I- HO2 -4 NO2 4- OH with minor contributions from (R24) 
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Figure 3.4: Sensitivities of OH to reaction steps predicted by WMR model. 


H + N0 2 -> OH + NO. 

Figure 3.5 shows the sensitivities of NO with respect to different steps. As revealed 
in the figure, NO is consumed mainly by its conversion to HN0 3 via the three-body 
reaction step (R51) OH + NO + M -4 HN0 2 4 - M. The next important reaction chain 
influencing NO consumption is step (R41) CO 4 - OH — ► C0 2 4 - H followed by step 
(R45) H 4 0 2 + M -4 H0 2 4- M. These two steps generate H0 2 which reacts with 
NO via step (R20) H0 2 + NO -4 N0 2 4 OH. Reaction (R43) 0 + 0 2 -^0 3 and the 
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Figure 3.5: Sensitivities of NO to reaction steps predicted by WMR model. 


further reaction of 0 3 with NO via (R19) NO + 0 3 NO 2 + 0 2 also contribute to 
NO reduction. Steps (RIO) OH + OH -t H 2 0 + O, (Rll) OH + H0 2 -» H 2 0 + 
0 2 , and (R30) HN0 2 + OH — > H 2 0 + N0 2 are the most important reactions for NO 
production. These reactions compete with the most dominant NO destruction step 
(R51) for OH. 

Figure 3.6 shows the computed sensitivity results of N0 2 . The most dominant 
reaction path for N0 2 production is through step (R41) CO + OH — > C0 2 + H, with 
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7. OH*0->H+O t 

10. 2OH*>H,O+0 
II.OH+HO.^HjO+O, 

11. N0+0 > «>N 0 , 4 - 0 , 

20. NO*HO,«»NO.+OH 
22. N0,*0«>N0+0, 

24. H*NO,«*OH+NO 

30. H NO, ♦OHk>H ,04-NO, 

41. CO+OHat>CO,+H 
43. 0*0,(4M)«>0,(+M) 

45. H*0,(*ll)*>H0 f (+M) 

43. 20H(+M)*>H,0,(4-M) 
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Figure 3.6: Sensitivities of NO 2 to reaction steps predicted by WMR model. 


subsequent reaction (R45) H + O 2 + M — > HO 2 + M and further reaction (R 20 ) NO 
+ HO 2 — » NO 2 + OH. NO 2 is consumed mainly by the three-body reaction (R53) 
N0 2 + OH + M HN0 3 + M. 

Figure 3.7 shows the sensitivity of ozone, O 3 , with respect to various reaction 
steps. O 3 is mainly produced by step (R43) 0 + 02 + M->0 3 + M and the 0 atom 
is supplied via reaction (RIO) OH + OH -4 H 2 0 + O. As OH is depleted mainly due 
to its reactions with NO, these two reactions become less important. O 3 is consumed 
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Figure 3.7: Sensitivities of O 3 to reaction steps predicted by WMR model. 


mainly by (R19) NO + 0 3 -4 N0 2 + 0 2 . As more 0 3 is entrained into the plume 
with time, reaction (R19) becomes increasingly important. 

Figure 3.8 shows the sensitivity of H 2 SO 4 with respect to different steps. The 
buildup of H 2 S0 4 is initiated by (R59) S0 2 + OH + M -4 H 2 S0 3 4 - M, and the 
product HS0 3 further reacts with 0 2 via fast reaction (R40) HS0 3 + 0 2 -4 H0 2 4 
S0 3 (not shown on the plot due to the small H 2 S0 4 sensitivity to this reaction). Fi- 
nally H 2 S0 4 is generated through (R39) S0 3 4 H 2 0 -4 H 2 S0 4 . Since OH is consumed 
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Figure 3.8: Sensitivities of H 2 S0 4 to reaction steps predicted by WMR model. 

mainly by NO and N0 2 via three-body reactions (R51) and (R53) and by self con- 
sumption reactions (RIO) and (R48) as depicted in Figure 3.4, the conversion of S0 2 
to H 2 S0 4 is expected to be low. This low conversion fraction will be demonstrated 
later in this section from the distribution of major sulfur species. 

Figure 3.9 to 3.11 depict the major species evolution along the plume axis. A 
second x-axis is also plotted on the top of these figures to show the corresponding 
plume age. The centerline evolution of major nitrogen species is shown in Figure 3.9. 
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At about 5 m past the nozzle exit plane, the turbulent stirring (macroscale mixing) 
becomes important and the corresponding plume age is about 3 ms. The primary 
exhaust products NO and NO2 stay nearly constant up to this point and then decrease 
slowly. In the early jet regime, small amounts of NO and NO2 are oxidized by OH 
(not discernible on this logarithmic scale) via three-body reactions (R 51 ) and (R 53 ). 
However, these reactions lead to significant amounts of nitrous acid, HNO2, and 
nitric acid, HNO3, as compared to their background abundances. A small fraction 
of HNO2 is consumed through (R 30 ) when the OH level is high. After OH becomes 
depleted (compare Figure 3 . 11 ), the chemical buildup of HNO2 and HNO3 ceases and 
these secondary exhaust products are diluted like tracers, similar to the fate of NO x 
in the late jet regime. At later stages of plume dispersion these acids may either 
become photolyzed under daytime conditions, or depleted by heterogeneous removal 
processes. 

In the early jet regime, the nitrate radicals, NO3, are generated via (R 52 ) NO2 -1- O 
+ M — > NO3 -I- M and (R 32 ) HN 0 3 -I- OH -> H 2 0 + NO3. The chemical production 
of NO3 ceases as soon as O and OH become depleted and its further evolution is 
governed by plume dilution. In the early jet phase, dinitrogen pentoxide, N2O5, is 
produced via (R 54 ) NO2 + NO3 + M — > N2O5 -I- M in contrast to the strong thermal 
decay reaction (R 55 ) N2O5 -t- M — > NO2 + NO3 + M. Hereafter the N2O5 abundance 
increases quickly by entrainment, and at the end of jet regime, N 2 0 5 already attains 
its ambient level. The NO3 level at the end of the jet regime is much lower compared 
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Figure 3.9: Centerline evolution of mole fractions of major nitrogen species. 

to its background abundance due to step (R54). 

In the the jet plume, the ratio of NO x (= NO + NO 2 ) to NO y (= NO* + HNO 2 
+ HNO 3 + 2 N 2 O 5 ), assumed to be unity at the nozzle exit plane, decreases slightly 
as NOy is built up. [NO x ]/[NO p ] is found to be approximately 0.975 at the end of jet 
regime. 

Figure 3.10 shows the centerline evolution of major sulfur species. The sul- 
fur reaction chain is initiated by oxidation of the exhaust S0 2 by OH via (R59) 



37 


S0 2 + OH + M -4 HS0 3 + M which immediately builds up HS0 3 . This is followed 
by fast reaction (R40) HS0 3 + 0 2 -> H0 2 + S0 3 , in which S0 3 is generated. Sul- 
furic acid, H 2 S0 4 , is then generated through (R39) S0 3 + H 2 0 -> H 2 S0 4 . Reaction 
(R59) is the rate-limiting step in this three-reaction sequence in producing H 2 S0 4 . 
Consequently, S0 2 oxidation is indirectly coupled with NO„, CO, and HO* reactions 

that compete for OH. The maximum conversion of S0 2 to H 2 S0 4 , ^ //2S °^ max j s 

(^so 2 )o 

less than 1%. After OH is depleted, the SO3 concentration reaches its maximum and 
then decreases due to plume dilution. H 2 S0 4 exhibits a similar evolution as S0 3 , 
except that it is shifted to slightly downstream distance. H 2 S0 4 losses due to binary 
H 2 S0 4 -H 2 0 nucleation is not considered at this time. It will be discussed in detail 
later in Chapter 4. 

Figure 3.11 depicts the evolution of major oxygen and hydrogen species along 
the plume axis versus downstream distance. Hydroxyl radical OH, being the most 
abundant primary HO* exhaust product, is responsible for the rapid production of 
other radicals. In the early jet regime, O is produced via (RIO) OH -I- OH — ► 0 
+ H 2 0, which further produces ozone through the three body reaction (R43) O + 
0 2 + M — ¥ O3 + M. In the late jet regime, mixing becomes important and the 
further increase of the O3 concentration is due to the entrainment of ambient ozone. 
However, as long as the NO level stays high, ozone is continuously consumed due to 
its reaction with NO via (R19) NO + 0 3 -¥ N0 2 + 0 2 and stays below its ambient 
level within the jet plume. The self-reaction (R48) OH + OH + M -> H 2 0 2 + M is 
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Figure 3.10: Centerline evolution of mole fractions of major sulfur species. 

the most important pathway for the production of hydrogen peroxide, H2O2, leading 
to a maximum abundance around the point of the strongest mixing. In the late jet 
regime, H2O2 decreases due to plume dilution. On longer timescales H2O2 may either 
act as a HOi sink through heterogeneous removal or as a slow source for regeneration 
of OH by photolysis. 

The buildup of hydrogen atoms H in the early jet regime is initiated by reaction 
(R41) CO + OH — > H + CO2. H further reacts with O2 via (R45) H + O2 + M — > 
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HO 2 + M which is the major production pathway of the hydroperoxy radical, HO 2 
Hydrogen molecules H 2 are subsequently produced via (R3) H + OH — ► O + H 2 and 
(R5) HO 2 + H — ► H 2 + O 2 , and by (R46) H + H + M-+H 2 + M when the H level 
is high. 

It is observed in Figure 3.11 that the OH concentration dramatically decreases in 
the intermediate stage of jet mixing between 1 and 10 m past the nozzle exit. As 
depicted from the OH sensitivity coefficients in Figure 3.4, the rapid removal of OH 
is mainly due to its reaction with NO and NO 2 via (R51) and (R53)and partially 
due to self reactions (RIO) and (R48). The reactions that sustain the O and H 
concentrations in the isothermal jet core are quickly reduced due to the depletion of 
OH radical. Hence the O and H atoms quickly decrease to negligible levels and only 
slowly recover in later stages of plume dispersion due to OH regeneration through 
photolysis. OH and HO 2 abundances sharply bend in the decaying jet flow field, 
which is due to the entrainment of HO 2 from the ambient air and its quick conversion 
to OH by reaction with NO through (R20) NO + HO 2 — > NO 2 + OH. 

Figure 3.12 shows the centerline temperature and major species distributions ob- 
tained with the detailed and the eight-step reduced mechanisms. It is assumed that 
the fluids in each computational grid are well-mixed. The centerline temperature 
stays nearly constant up to the point where the ambient air begins to mix with the 
centerline gaseous species. Temperature decreases rapidly and reaches near ambient 
level at the end of jet plume. Predictions of the temperature and major species evo- 
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Figure 3.11: Centerline evolution of mole fractions of major hydrogen and oxygen 
species. 

lutions by the eight-step reduced mechanism is shown to be in excellent agreement 
with those by the detailed mechanism. 

Figure 3.13 shows the results obtained under the same conditions as in Figure 3.12 
but with the subgrid mixing accounted for by the modified Curl’s mixing model. 
When the interaction of subgrid mixing with chemistry is considered, the eight-step 
reduced mechanism still gives results that are in good agreement with those predicted 
by the detailed mechanism. Simulation with the reduced mechanism is about 3.4 times 
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faster in computer CPU time compared to the computation with the detailed one for 
the near-field jet plume. 

Figure 3.14 compares the centerline major species distributions with and without 
the interaction of subgrid turbulent mixing and chemistry. The differences between 
these two profiles are relatively minor. The predicted results with the subgrid mixing 
modeled by the modified Curl’s mixing model show slightly lower levels of species 
concentration compared to the results with the subgrid well-mixed. 

There have been no measurements of major species concentrations carried out in 
the jet plume regime. Fahey et al. [1995] carried the in-situ measurement with an 
ER 2 research aircraft flying in the lower stratosphere during the NASA Stratospheric 
Photochemistry, Aerosol and Dynamics Expedition (SPADE) program. To compare 
our model predictions with their measured values, we perform calculations for an ER 
2 aircraft with flight conditions and initial exhaust properties given by Anderson et 
al.[1996]. Figure 3.15 compares the ratios of various chemical species concentrations 
above the ambient background levels for an ER 2 aircraft from the present predic- 
tions at 1 km to UNIWAKE results at 20.2 km by Anderson et al. [1996] and to 
the measured values at 126.0 km by Fahey et al. [1995]. The present calculations 
present results at a much smaller distance downstream compared to the UNIWAKE 
results and the SPADE data. However, since chemical kinetics is negligible after 
100 m downstream and the evolution of species concentration is due almost entirely 
to mixing, the ratios of AN0 2 /AN0, ANO x /ANO y , and ANO y /AC0 2 are gen- 
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erally independent of the dilution level and thus tend to stay constant downstream. 
The comparison shows good agreement among different predictions and the SPADE 


measurements. 
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Figure 3.12: Centerline temperature and major species distributions with the detailed 
(lines with squares) and the reduced (lines with reversed triangles) mechanisms. Sub- 
grid is assumed to be well-mixed. 
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Figure 3.13: Centerline temperature and major species distributions with the detailed 
(lines with squares) and the reduced (lines with reversed triangles) mechanisms. Sub- 
grid mixing is modeled by modified Curl’s mixing model. 
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Figure 3.14: Centerline major species distribution by pdf jet simulation, with different 
treatment of subgrid mixing model. Lines with squares are results with subgrid mixing 
well-mixed, and lines with reversed triangles are results with subgrid mixing modeled 
by modified Curl’s mixing model. 
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Figure 3.15: Comparison of the ratios of various chemical species concentrations over 
the ambient background levels from the present predictions at 1 km to UNIWAKE 
results at 20.2 km and to measured values at 126.0 km. 
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Chapter 4 

Aerosol Dynamics 

4.1 Introduction 

The emission and formation of new aerosols caused by subsonic aircrafts and their 
subsequent evolution in the upper troposphere and lower stratosphere have become 
an important research issue. Global impact of fleets of aircraft depend on how much 
H 2 SO 4 -H 2 O aerosols are generated in the near-field aircraft exhaust [Weisenstein, 
1996]. In this regime, not only major gas- phase chemical reactions take place, but 
also a large number of small-sized volatile H2SO4-H2O aerosols are generated due to 
binary homogeneous nucleation. Soot activation as water condensation nuclei also 
occurs in this regime. Preliminary analysis of data collected under the SUCCESS 
(SUbsonic Aircraft: Contrail and Cloud Effects Special Studies) program {NASA, 
1996] has indicated a high level of small volatile aerosols whose number density varies 
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with the mass loading of fuel sulfur [Ballenthin et al., 1996]. 

In-flight field measurements of aerosols by the instrumentation aboard the ER2 
aircraft in plume encounters of the Concorde wake [Fahey et al., 1995] have reported 
between 12% and 45% conversion of fuel sulfur to sulfuric acid (assuming the diameter 
of sulfuric acid aerosols is 10 nm) after approximately 16 minutes at 16.2 km from the 
nozzle exit plane. Ground-based measurement indicates sulfur oxidation greater than 
0.4% at the early plume stage [Frenzel and Arnold, 1994]. Modeling studies of aircraft 
plume chemistry have predicted that the conversion of fuel sulfur to sulfuric acid is 
limited to 0.5%-2% [Karcher et al., 1996a, Brown et al., 1996a, Wang and Chen, 
1997]. These numerical estimates are based on the assumption that fuel sulfur is 
emitted from the engine as SO 2 which is subsequently converted to sulfuric acid by gas 
phase reactions in the exhaust plume. The discrepancy seen between the experimental 
estimate by Fahey et al. [1995] and the numerical results raises concern on the validity 
of the above assumption. Hence, field measurements of sulfur speciation at the nozzle 
exit plane are suggested to resolve the disagreement. 

Using numerical modeling, Brown et al. [1996b] have recently estimated the ex- 
tent of fuel sulfur converted to volatile aerosol precursors during combustion and the 
subsequent turbine flow for the subsonic ATTAS and the supersonic Concord air- 
craft engines. The results of their analysis suggested that the conversion rate of fuel 
sulfur to S(VI)(S0 3 and H 2 S0 4 ) is between 2% and 10% in the engine. The large 
uncertainty in the estimated values is due to the limitations of the models used in 
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that study to characterize the oxidative pyrolysis of typical jet engine fuels and the 
three-dimensional fluid dynamics for real combustors and turbines. The lack of high 
temperature kinetic data for important sulfur oxidation reactions also contributes to 
the uncertainty. 

The laboratory studies by Wyslouzil et al. [1994] with samples of carbon black 
have indicated that carbon is hydrophobic. This result suggests that soot particu- 
lates emitted in aircraft exhausts may be poor condensation nuclei in the absence of 
activating precursors. Both Wyslouzil et al. [1994] and Rogaski et al. [1996] have 
observed in their experiments a substantial increase in water uptake on carbon in the 
presence of sulfuric acid. Therefore, the potential role of sulfuric acid formed through 
sulfur oxidation is currently of particular interest to soot activation. 

In this chapter, we will consider the microphysical processes including binary ho- 
mogeneous nucleation, aerosol coagulation, adsorption, and heteromolecular conden- 
sation. These processes are coupled with gas phase chemical reactions and turbulent 
mixing in an axisymmetric flow. Compared to the modeling work of Brown et al.’s 
[1996a; 1996c], the present work assesses the impact of SO 3 emissions on the formation 
and evolution of aerosols and studies the condensation processes with a kinetic model. 
Karcher [1997] studied the near-field microphysical processes with a box model, in 
which the aerosol processes were not coupled to the flow field and no spatial resolution 
was offered. The mixing and cooling rates along the plume trajectory were chosen 
to represent the evolution of the bulk of the exhaust. Currently, there are only few 
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aerosol modules implemented in global chemistry models [Bekki and Pyle, 1992; Tie et 
al., 1994; Weisenstein et al., 1997]. However, the formulations of aerosol microphysics 
used in global models are too simplified to be applied to near-field aircraft-related 
aerosol simulations. 


4.2 Aerosol Model Structure 

The time evolution of number density due to homogeneous nucleation, coagulation, 
and condensation can be described by the following equations 


j _ /dN V j dN V j. 

dt — J * lom ' fit Jcoag * \ ^ Jcond 

(4.1) 

dN c , . dN Cii v . dN Cyi . 

dt dt )coa9 ( dt )amd ' 

(4.2) 


where Ni is particle number density representing particle number in size bin i per 
cm 3 of air with the unit of #/cm 3 , and J denotes the nucleation rate. The subscript 
v indicates volatile H 2 SO 4 -H 2 O aerosols, and the subscript c indicates either initially 
emitted bare soot particles or later evolved soot-cored droplets. It is computationally 
efficient to combine the number density of bare soot and that of soot coated with 
H 2 SO 4 -H 2 O as one variable, which will be justified in section 4.5. Given the number 
density N, the mixing ratio x can be computed by x = N/N a , where N a is the number 
density of air ( N a = P/{k B T) where k B is Boltzmann’s constant). Details of each 
term in equations (4.1) and (4.2) will be described in the next section. 
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4.3 Aerosol Model Formulation 

The size range of particles is discretized into a number of j max radial size bins. 
Each set of size bins is based on a geometrically increasing particle volume. The bins 
are centered at volume 


v 3 = v rat v 3 _ x 


(4.3) 


with Vrat > 2. Correspondingly, the particle radius r 3 is located at 

r j = (VRAT) l/3 rj-\- 

The bin widths increase according to 

Vrat - 1 


AV 3 = 2 Vj 


V RAT + 1 


/ o \ 1 / 3 

- (vwtt) (v - - 


The upper and lower limits of each particle size bin in terms of radius are 

( rat \ 1/3 

rj '" = (vwTTj r ' 

= (vwTTj r * 

The total number of bins covering the particle size range from r\ to r 3 is 

_ lnfo/rQ 8 

Jmax — i "r 

m Vrat 


(4.4) 


(4.5) 

(46) 


(4.7) 

(4.8) 


(4.9) 


In this study, we use ri = 0.3 nm which is about the size of the sulfuric acid embryo, 


V R at = 2, and j max = 45 which leads to r 3max ss 8 fim by doubling the volume for 


consecutive bins. 
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The typical size distributions of nascent primary soot particles are assumed to be 
log-normal, with median radii of 15-30 nm, modal widths of 1.4-1 .6, and number 
densities of 10 6 -10 7 cm" 3 [Hagen et al., 1992; Petzold and Schroder, 1998]. The 
log-normal distribution for the size of soot particles can be expressed as 


n(lnr) = 


N t 


exp 


(lnr — lnr 9 ) 2 
2 In 2 a„ 


(4.10) 


\/27r In o g 

where n(lnr)dlnr is the number of particles per cm 3 of air in the size range from lnr 
to In r -I- din r, N t is the total number of particles of all sizes per cm 3 of air, f 9 is the 
median radius, and o g is the geometric standard deviation. The number density of 


particles in the jth size bin centered around the radius is given by 


dN = n(lnr)dlnr — n ^ n r ) ( j r 

r 

N t 


y/2nr In a g 

4.3.1 Homogeneous Nucleation 


exp 


(lnr — lnr 9 ) 2 
2 In 2 o a 


(4.11) 


dr. 


The classical nucleation theory is based on the equilibrium state of critical-sized 
liquid droplets, or embryos, in contact with their parent vapors. The theory gives the 
nucleation rate in the form of 


J = A exp (—AG*/(R U T )) , (4.12) 

where J is the nucleation rate with the unit of embryos/ (cm 3 -sec), A is the kinetic 
prefactor, AG* is the change of the Gibbs free energy during the phase change for 
the formation of an embryo, and R^, is the universal constant. The exponential 
dependence of J on AG* implies that AG* must be accurately calculated. 
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The change of the Gibbs free energy for the formation of a liquid embryo in contact 
with a binary vapor mixture can be written as [Flood, 1934; Reiss, 1950; Doyle, 1961] 

AG = n - n\) + n 2 (/z 2 ~ f4) + 4tt r* 2 a, (4.13) 


where rii is the number of moles of component i in the embryo, n\ is the chemical 
potential of the ith component in the liquid solution with the unit of J/mole, pf is 
the chemical potential of the ith component in the gas phase, a is the surface tension 
with the unit of J/cm 2 , and r* is the radius of an embryo. The first two terms on the 
right side represent the energy that is gained by forming the liquid phase and the last 
term is the energy required to create the surface of the liquid droplet. The change of 
chemical potential in transforming from gas to liquid phase for each component is 

= In 


fil ) 
p^j ' 


(4.14) 


where Pi represents the partial pressure of gas-phase species i and P°^ ot is the equi- 
librium vapor pressure of species i over a flat surface of the mixed solution. 

The critical free energy AG* is determined as the value of AG at the saddle point 
in the two dimensional energy surface. The conventional way is to solve iteratively 
a pair of generalized Kelvin equations [Mirabel and Katz, 1974] representing the 
derivatives of AG with respect to ni and n 2 : 



(4.15) 


However, this method involves tedious calculations in a two-dimensional parameter 
space. Recently, Zhao and Turco [1995] have devised an efficient method for locating 
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the saddle point (critical Gibbs free energy) that reduces the search for the roots of 
the generalized Kelvin equations in a two-dimensional space (ni,n 2 ) to a search in 
two one-dimensional spaces. They transformed the coordinate system from (ni,n 2 ) 
to (u, V), that is, the weight fraction of sulfuric acid and the volume of a solution 
droplet. They are defined as 


M 2 n 2 

Mini + M 2 n 2 


(4.16) 


and 


V = 


Mini + M 2 n 2 


(4.17) 


where Mi is the molar weight of species i and p is the density of the solution. Then 
a revised version of the generalized Kelvin equations is derived from the derivatives 
of AG with respect to V and u> for a binary system as 


dV_o_ 
lA~ Y' 


(4.18) 


A/i 2 


dlnp 

1 + -—U 

du 


M2 

Mi 


Api 




= 0, 


(4.19) 


where 


Y = -p 


u> . 1 — uj . 

A + "mT A, “ 


(4.20) 


and A = 4nr 2 is the surface area of an aerosol. 


The advantage of the transformation is that equation (4.19) is independent of 
radius so that u* can be obtained iteratively by solving this equation alone. Equa- 
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tion (4.18) can be simplified as 


r* = 


2a (lo*) 

V>')' 


(4.21) 


Hence, once u* is determined, r* can be obtained immediately with equation (4.21). 
Using equation (4.13), the Gibbs free energy at the saddle point is 


AGw = -7ra(cj*)r* 2 . 


(4.22) 


The homogeneous nucleation rate with the unit of embryos/(cm 3 -sec) is 

■ «•“> 

where /? 2 is the gas kinetic impingement rate of H2SO4 molecules towards the embryo 
surface, 0 2 = N 2 (kBT/27rm 2 ) l/2 with the unit of molecules/(cm 2 -sec), and N\ is the 
number density of water molecules in the gas phase. The prefactor, A = 47rr* 2 /? 2 7Vi, 
is valid for the usual situation when N\ » N 2 . 

Finally, the critical numbers of H 2 S0 4 and H 2 0 molecules in an embryo are then 
computed from 


Jhom = 47T r 


* 2 02 Ni exp ^ 


N; = p(u\T)\r(l-w’)/mi 

(4.24) 

N; = p(u-,T)Vu-/m 2l 

(4.25) 


where the volume of the embryo is V* = 4nr* 3 /3, mi and m 2 are the molecular masses 
of H 2 0 and H 2 S0 4 , and p is the mass density of the liquid solution droplet. 
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4.3.2 Coagulation 


Coagulation model is based on Toon’s [1988] approach which conserves particle 
number and volume. The change of number density of volatile H2SO4-H2O particles 
due to coagulation can be described by 


(^U, = + N.t £(1 - QlCvN.j 

j<i * 

+N v ,i - 1 El 1 " - N «-< E 


j<*-l 




-JV.,, E K V N ‘J < 4 ' 26 ) 

J 

where Kij — K }1 is the Brownian coagulation coefficient between particles with size 
indices i and j. The subscript v denotes the volatile H2SO4-H2O particles and c 
denotes the soot or the soot coated with H2SO4-H2O. 

Because bin i has twice the volume of bin i — 1, collisions between bin i and bin j , 
where j is less than i, result in particles of volume intermediate between bin i and bin 
i + 1 . The volume ratio Vj/V, is the fraction added to the larger bin by this collision. 

Boundary conditions in coagulation are imposed to conserve mass. Instead of 
allowing particles to leave the largest size bin due to coagulation, those particles are 
kept in the largest size bin. Special conditions also apply to the production terms in 
the first size bin. It is assumed that there is no particles of size smaller than the first 
size bin so that coagulation cannot produce particle of size j — 1. 

The Brownian diffusion kernel is calculated using the interpolation formula by 



57 


Fuchs [1964] 




47r(r, + rj)(Di + Dj) 


rj + rj 


4(A + Dj) 


r, + Tj + (<5, 2 + 5 2 )V 2 - (V 2 . + + r .) 


1. 

/? 


(4.27) 


where r* and r ; are the radii of particles in the ith and jth size bins and /? is the 
sticking coefficient which is assumed to be unity in the present work. D x and Dj are 
the particle diffusion coefficient that can be calculated from the Einstein relation 


Di = ksBjT, 


(4.28) 


where ks is the Boltzmann constant. T is the absolute temperature and B is the 
mobility given by 


B,= 


1 

6nr]r t 


1 + 1. 246/^71* + 0.42 exp 



(4.29) 


where rj is the viscosity of air and Kn is the Knudsen number defined as the ratio of 
mean free path of air and particle radius 


Kn x = — . 

Ti 


(4.30) 


v Pi i is the mean thermal velocity of a particle 


/8 k B T\' n 


( 4 . 31 ) 


where rrii is the mass of the particle. <5, is a correction factor given by 

, (2r, + A*) 3 - (4rf + „ 

». = t— : 2ri, 

6?": Apj 


(4 32) 
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with the mean free path of a particle in size bin i being 

(4.33) 

It can be shown that the coagulation coefficient increases with the difference be- 
tween rj and r 2 and it reaches minimum for coagulation between particles with the 
same size. Coagulation between the soot particles is not important due to their 
low amount compared to the large amount of H 2 S0 4 -H 2 0 particles, and therefore is 
neglected in the wake-regime studies. This is because the time scale required to intro- 
duce a significant change in the distribution function for these aerosols is longer than 
the time available in the near-field plume. The characteristic time for coagulation is 
estimated as [Flagan and Seinfeld, 1988] 

tcoas = ^H' 

where n is the number density of particles and the typical value for the coagulation 
constant is H = 10 -9 cm 3 /sec. This implies that the number density must lie between 
10® and 10 10 cm -3 in order for coagulation to be important in the near-field exhaust 
flow. The number density of H 2 S0 4 -H 2 0 embryos formed through binary homoge- 
neous nucleation is sufficiently large for embryo-embryo and embryo-soot interactions 
to be important on a timescale of a few seconds. However, the soot particulates 
emitted in the exhaust flow are typically present at much lower levels. 

Scavenging of H 2 S0 4 -H 2 0 aerosols by the soot particles is an important soot 
surface chemical activation pathway. Wetting the soot particles by this mechanism 
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is assumed to have no effect on the size of the soot because it can only lead to a 
monolayer coverage. However, the heterogeneous condensation on the liquid coating 
of soot is coupled to radial growth and might cause a shift in the size distribution. 

4.3.3 Soot Activation 

Soot particles can be activated as water condensation nuclei by adsorption of 
oxidized sulfur species (H2SO4 and SO3) and by scavenging (coagulation) of volatile 
sulfate aerosols. The fraction of soot surface that is covered by oxidized sulfur species 
is defined as 9{ = >1(503 + where >1 S|< is the surface area density of soot 

in size bin i, which is A-nr]^N cd . The surface coverage due to the above two activation 
pathways is 9 , = 0 ads>i + 9 sca>x . 

Karcher [1998] presented a kinetic formulation which gives the maximum sulfur 
adsorption on soot. Since r sd < 100 nm, vapor adsorption on fresh soot takes place 
in the gas kinetic regime. Therefore, the corresponding rate of change of 9 is given 
by 

a 1* ,t = o:-(Cso 3 + C Hi so<) x 6 02 x 10 23 (4.34) 

ut 4 CTq 

where a is the sticking probability which is set to zero for T > 420 K and equal 
to unity for T < 420 K, v is the thermal speed for the SO3 and H2SO4 molecules 
evaluated at their reduced mass, C is the molar density of gas phase species, and 
ffo«5x 10 14 cm -2 is the approximate number of active sites available on the soot 
surface. This kinetic formulation gives the maximum adsorbed sulfur mass and it 
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represents the upper bound to the binary H2SO4-H2O heterogeneous nucleation rate 
on the soot surface. 

The change of surface coverage for a soot particle in size bin i by scavenging 
volatile droplets is determined by 


de„ 


Jmoz 


dt 


' a = Y,k( 

3-1 


r»,ur vJ )N vJ nr:,; x 




i-», 


( 4 . 35 ) 


where K is the coagulation coefficient and N v is the number density of volatile droplet. 
The gas phase losses of SO3 and H2SO4 due to adsorption is oc (1 — 9 ), i.e., they are 
proportional to the available “dry” soot surface area 


(^r) = -o \c m 4<,JV C , ( (1 - e,), (4.36) 

V / ads \ 

where m denotes H2SO4 or SO3. This model also accounts for the condensation of 
H2SO4 and H 2 0 on the part of the soot surface that has already acquired the liquid 
coating. Consequently, condensation rates on soot and related gas phase losses scale 
with 6 . 


4.3.4 Condensation 

Condensation can occur on the liquid volatile aerosols and on the partial liquid 
coating on the soot surface. The equilibrium vapor pressure of H 2 S 0 4 over the solution 
is so small that once H2SO4 molecules are condensed they can not re-evaporate. In 
contrast, droplets may take up or evaporate H2O, depending of the evolution of 
temperature and the partial pressure of H 2 0 in the plume. 
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The condensation formulation follows the derivation by Fukuta and Walter [1970], 
which was based on the idea of Fuchs [1959]. The space outside of a particle is divided 
into two zones, the layer r < R < r+6 and the region beyond R > r+<5 where 6 « A 0 , 
the molecular mean free path of air. The conventional diffusion equation, which is 
based on the assumption that vapor density in air is continuous, is valid for R > r + 8, 
while gas kinetic theory is applied within the layer. The condition of continuity of 
vapor flux across the surface R = r + S is invoked to complete the description. 
Consequently, the change in moles of gaseous species m due to condensation on or 
evaporation from an aerosol with radius rj is 

t (Pm Pm,i,sat\ i. n ~,\ 

— = { — hT ~ ) , (4.37) 

where the modified diffusivity is 



Dy 

r\ 4A, 

+ S r t av m 


(4.38) 


with q being the mass accomodation coefficient. For condensation onto liquid surfaces 
considered in the present model, a is taken as unity. The mean thermal velocity of 
gaseous species m is given by 


^771 



(4.39) 


Without the diffusivity modification, equation (4.37) is in the conventional Maxwell’s 
form, which would strongly overestimate the growth rate for submicron particles. For 
soot particles, since condensation can only take place on the activated fraction 6 of 
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the surface, the right hand side of equation (4.37) needs to be multiplied by 6 when 
calculating condensation on soot particles. 

The saturation pressure is obtained by the Kelvin equation 


p 

1 m,t,, 


sat 


m>oo } sat 



(4.40) 


where P m ,i )SO t is the saturation vapor pressure of gaseous species m over a spherical 
droplet of radius r<, and P m ,oo,jot is the saturation vapor pressure of species m over a 
flat surface of a liquid with the same composition as the droplet. The Kelvin equation 
(4.40) states that the saturation vapor pressure over the surface of a droplet is greater 
than that over a flat surface. The smaller the droplet, the larger is the saturation 
vapor pressure. Consequently smaller droplets tend to evaporate, while larger ones 
tend to grow. The effect of curvature is important for droplets less than 0.1 ^m. 

Denoting the total radial growth rate of a particle (volatile droplet or soot coated 
with H2SO4-H2O) in size bin j due to condensation by dr 3 /dt, the flux of particles 
growing from tj to ry + A r } per unit time and per cubic centimeter of air is 

r .= *L±i 
J A Tj dt ’ 



dv 

where the expression of —■ will be given later in this section. The rate of change 

dt 


of particle abundances due to condensation ( Gj > 0) and evaporation (Gj < 0) will 
not affect the total number of each particle type. Therefore, the net rate of change of 
each particle type in size bin j over a time step can be expressed as [Karcher, 1998] 


A Gj = max{0, Gj- 1 } + | min{G J+ i,0}| - |Gj|. 


(4.42) 
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The first term on the right hand side represents the possible gain of particles in size 
bin j from size bin j — 1 upon condensation. The second term denotes the gain of 
particles from size bin j + 1 due to evaporation. The last term is the loss of particles 
in size bin j by growth or evaporation. To conserve the total particle numbers in each 
type, boundary conditions at size bin j = 1 and j = j max are set so that particles are 
kept within the size range [r m j n ,r moi ]. 

Using equation (4.42), the term ( dN/dt ) amd in equations (4.1) and (4.2) can be 
expressed as 


/ dNj\ _ AGj __ AGj 

V dt N a - ( P a /k B TY 


(4.43) 


The change of radius due to condensation or evaporation is 


dr 

dt 


l _ XT ( dr, \ 

^ m \d n mj / 


dn. 


m,j 


dt ’ 


m = H 2 0,H2S0i, 


(4.44) 


where the partial derivative can be obtained from mass conservation in the liquid 

drj M m 

sr 3 - = 1 — ^s- ( 4 - 45 ) 

on m j 4 tt PsoiTj 

The condensation and evaporation processes control the weight fraction of H2SO4 

in the liquid solution, Wh 2 so 4 > on which the thermodynamic quantities depend. This 

especially holds for the saturation vapor pressure which in turn determines the con- 
dfi 

densation rate —— . Therefore it is necessary to determine W HiSOi ■ FYom the 
definition 


U'// 2 so 4 j 


MHiSoJVjhso* 


m = H 2 SO 4 , H 2 0 , 


(4.46) 
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where 


M m N mj = 4np sol r*/3, (4.47) 


and taking derivative of equation (4.46) with respect to time and combining equations 
(4.44) and (4.45), we get 

dWj _ 3 

dt rj 


f-iu 

\ dn m j ) 


d n mj _ drj 
dt J dt 


m = H 2 O, H 2 S0 4 . 


(4.48) 


4.3.5 Summary of Coupling between Chemistry and Aerosol 
Dynamics 


The evolution equations of gas-phase H 2 SO 4 and H 2 O will have additional sources 
or sinks due to the microphysical processes mentioned in the previous sections. The 
rate equations for H 2 O and H 2 SO 4 are modified as follows 


^ =*1 - - E - E - E Q i c ™X > r « w «< 1 - *), 


i=i 


(4.49) 


tmoz 


^ « 2 - - E - E “;On E - «,). 


»=i 


(4.50) 


where the subscript 1 denotes H 2 O, and 2 denotes H 2 SO 4 , C is the species gas-phase 
concentration with unit of mole/cm 3 , and u> is the gas-phase chemical reaction rate. 
Jhom is the homogeneous nucleation rate given by equation (4.23), n* m is the moles of 



gaseous species m in a H 2 SO 4 -H 2 O aerosol of critical size, and N, is number density 
of particles in size bin i with the subscript v as volatile H2SO4-H2O aerosol and 


dn^ 

dt 


the subscript c as soot coated with H 2 SO 4 -H 2 O. Terms due to condensation, 
dn 2 i 

and ’ , can be obtained from equation (4.37). The last term accounts for soot 

dt 

activation. Equations (4.49) and (4.50) reveal the coupling of microphysical processes 
with chemistry kinetics through concentrations of H 2 SO 4 and H 2 0. Microphysical 
processes also couple with turbulent flow field from mixing of the exhaust plume with 
the ambient air. 


4.3.6 Thermodynamic Properties 

Various thermodynamic quantities are needed for simulation of the aerosol dy- 
namics. Described in this section are details of methods used to calculate or estimate 
these properties. 


Mean free path: 


At the standard condition, Pq = 1.013 x 10 5 Pa and To = 293.15 K, the molecular 

k^T T 

mean free path of air is A„o = 6.6 x 10 -6 cm. Since A 0 = — = oc — , for other 

V2PncP m P' 


levels of P and T, 
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Diffusivity: 


The diffusivity of water vapor D v in air is obtained from the extrapolation by 
Pruppacher and Klett [1997] with the following best-estimated relation for tempera- 
tures between -40°C and 40°C 



(4.52) 


with T 0 = 273.15 K, P 0 = 1013.25 mbar, and D v in cm 2 /sec. As D v a v\ a and 
v oc l/VW v , we get D v oc 1 /\/W v . Therefore, D Vt n 3 so 4 can be obtained through 



(4.53) 


Surface tension: 


The surface tension of H2SO4-H2O solutions has been reported by Sabinina and 
Terpugow [1935] for 15 levels of H2SO4 content. Data were reported for temperatures 
between 10°C and 50°C in 10°C intervals. The reported data are fitted with an 
eighth-order polynomial at each temperature. The surface tension at temperatures 
between tabulated values is then determined by the linear interpolation. 


Solution density: 

The solution density is obtained from tabulated values [Perry, 1963]. Tabulated 
data were given for temperatures between 0°C and 100°C and for H 2 S0 4 weight 
percent between 0 and 100. The tabulated data at each temperature are fitted with 



an eighth-order polynomial. The density at temperatures between tabulated data is 
also computed by the linear interpolation. 


Equilibrium partial vapor pressure: 


The equilibrium partial vapor pressure over the binary solution H2SO4-H2O is 
presented by Taleb et al. [1996] in a very simple form, which is valid over the entire 
range of concentrations for temperatures in the range 190-298 K. In terms of a Van 
Laar type equation, the activity coefficients 7! and 72 of components 1 and 2 in a 
binary system can be expressed as [Taleb et al., 1996] 


T log 7! = 
T log 72 = 


A x X\ 


[X 2 + B.Xtf 
A 2 Xl 


(4.54) 


[Xi -(- B 2 X 2 ] 2 ’ 

where T is the temperature in degree K, X is the mole fraction, and A t and B x are 
coefficients for a given binary solution. Based on the data of Zeleznik [1991], the 
following expressions are used in equation (4.54) 


A x = 2.989 x 10 3 - 


+ 


. 2.147 x 10 6 2.33 x 10 8 

r 

,4 2 = 5.672 x 10 3 - 


ji 2 


4.074 x 10 6 4.421 x 10 8 

^ 1 “ 


J>2 


(4.55) 


Bi = ^- = 0.527, 

n 2 


where A x and B x are for water, and A 2 and B 2 are for sulfuric acid. 


The definitions of the activity coefficient 7, and the activity a, are given by 
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7t = 


a L 
Xi 


Psol,i 

XJ\i' 


(4.57) 


where P go i,i is the equilibrium partial pressure of component i over a solution contain- 
ing i and Po t i is the equilibrium pressure of pure i. It is straightforward to compute 
Psoi,i provided that P 0 ,» is known, which is generally the case, through 


Pi iol,i — i T'o.i - 


(4.58) 


The equilibrium vapor pressure of pure water expressed in millibars for the temper- 
ature range between 220 and 320 K is taken from the compilation of Tabata [1973] 

, „ . . 1827.17834 71208.271 

log P 0 ,h 2 o (mbars) = 8.42926609 - ^ • (4.59) 

And the equilibrium vapor pressure of the pure sulfuric acid is given by Ayers et al. 
[1980] 

lnP 0 ,// 2 so4(atm) = -10156/T + 16.259. (4.60) 


Mixture molar volume: 

The molar volume of a two-component solution is given by 

M (1 — X2)M\ + X 2 A /2 
P P 

where xi is the mole fraction of component 2. For a two-component mixture, the 
partial molar volume of component 1 or 2 in a solution can be obtained through the 
tangent line at a given mixture composition on the composition curve v versus X 2 (can 
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be obtained experimentally) at the intercept x 2 = 0 or x 2 = 1 as shown in Figure 4.1. 
Mathematically, the partial molar volumes are expressed as 

dv 

V\ = v — — — x 2 (4-62) 

dx 2 


and 


v 2 = v — 


dv 

dxo 


(*2 - 1 ), 


such that (1 — x 2 )v\ + x 2 v 2 = v. 


(4.63) 


4.4 Threshold Condition for Contrail Formation 
from Aircraft Exhaust 

Contrails (condensation trails) are observed behind an aircraft only under certain 
conditions depending on various factors. The threshold criterion for contrail formation 
was first developed by Appleman [1953]. He explained that contrails are expected to 
form in the engine exhaust plumes when isobaric mixing between the hot and humid 
exhaust gases and cold ambient air results in a mixture reaching saturation with 
respect to water. However, Busen and Schumann [1995] observed contrails at ambient 
temperatures well above the Appleman’s threshold temperature. They explained this 
discrepancy by considering a fraction of combustion heat transformed to the trailing 
vortex system in the form of kinetic energy. They identified this fraction with the 
overall propulsion efficiency of the aircraft, i.e., the fraction of combustion heat that 
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is converted into work to propel the aircraft 

V = FV/ (Qm F ), (4.64) 

where F is the thrust, V is the airspeed, Q is the specific combustion heat, and rkp is 
the fuel flow rate. Hence, the exhaust gases contain only (1 — r])Q of heat per mass of 
fuel instead of all the heat assumed by Appleman [1953]. Consequently, the threshold 
temperature for contrail formation is higher than that if the exhaust gases are assumed 
to receive all the combustion heat directly. The overall propulsion efficiencies, rj, 
for modern bypass turbofan engines is usually higher than those for older turbojet 
engines or turbofan engines with low bypass ratio. The kinetic energy of the jet plume 
is shown by Schumann [1996] to have very little impact on the threshold temperature 
although the kinetic energy of the trailing vortex system induced by the aircraft is 
important. Hence it is justified to assume a stagnant plume when calculating the 
threshold criterion for contrail formation. 

Figure 4.2 illustrates schematically how the threshold ambient temperature can be 
determined in a partial water pressure versus temperature diagram. Saturated states 
corresponding to liquid water and ice are denoted by P, a t,uq and P iat ice , respectively. 
Mixing of the aircraft exhaust with ambient air follows along the line connecting the 
state of the exhaust gas at engine exit and that in the ambient air. The coordinate 
corresponding to the condition at the engine nozzle exit is (T exit ,x/i }0< exitPa) where P a 
is the ambient total pressure, which is on the far right upper corner outside the range 
plotted. The ambient condition corresponds to the point (T a , RH a • Psat{T a )) where 
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RH stands for relative humidity and the subsript a denotes the ambient condition. 
Visible contrails form when the the plume reaches supersaturation with respect to 
water [Jenson et al., 1998]. In Figure 4.2, this corresponds to the conditions when the 
mixing line crosses the saturation curve with respect to liquid water, P sa t,uq (given 
by equation (4.59)). Thus the threshold condition is defined as the situation when 
the mixing line is tangent to the liquid saturation curve. Also plotted in the figure 
is the saturation vapor pressure with respect to ice, P 3a t,ice(T ) which is given by 
InPiat.ice = 24.313 - 6146. S/T [Tabazadeh et al., 1997]. 

The gradient of the mixing line is approximated by [Schumann, 1996] 


_ ElH 2 OCpP a 

eQ(l~v)' 


(4.65) 


where e = MH 2 o/M a 
calculated through 


= 0.622. The temperature at the tangent point Tim can be 


dP 3a t(TiM) _ ^ 
dT = G 


(4.66) 


Equation (4.66) is nonlinear and it can be solved for Tim through Newton’s iteration. 
A good guess to initiate the iteration is given by [Schumann, 1996] 


T lm = -46.46 + 9.43 ln(G - 0.053) + 0.720[ln(G - 0.053)] 2 + 273.15, (4.67) 


where T is in unit K and G in is unit Pa/K. From the definition of the ambient 
relative humidity RH a and that of the gradient of the mixing line, the threshold 
ambient temperature for contrail formation is determined by 


TlC — Tlm — [ Psat{TLM ) — PHa ■ Psat{Tic)]/G. 


(4.68) 
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Equation (4.68) is also an implicit equation which requires Newton’s iteration. Tlm 
is a good initial guess to start the iteration. Contrails are expected to form when the 
ambient temperature is less than Tlc > for example, the temperature at point A as 
illustrated in Figure 4.2. The threshold ambient temperature can also be determined 
by the ice saturation state as denoted by T jc- As the specific volume increases when 
liquid water freezes, the saturation pressure of ice is higher than that of liquid water; 
consequently, Tjc > Tlc- 

In-situ measurements of the ambient temperature and humidity made by Schu- 
mann et al. [1996] and by Jensen et al. [1998] during the SUCCESS experiment 
suggest that liquid saturation in the plume is required for visible contrail formation. 
The fact that the ambient temperature needs to be below Tlc may indicate that liq- 
uid water droplets must be activated and subsequently freeze in order to generate a 
sufficient number of ice crystals for a visible contrail. According to the above in-situ 
observations, visible contrails do not form when Tlc <T < Tjc- In this temperature 
range, the number of ice crystals nucleated may be too low to form ice contrails that 
are visible. 

In this thesis work, the formation of ice particles by freezing and their subse- 
quent evolution are not modeled. This is justified for situations when contrails are 
not formed, i.e., when ambient temperature is above the threshold ambient tempera- 
ture. The present analysis also sets a stage for further research on ice particles when 


contrails are formed. 
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4.5 Results and Discussion 

4.5.1 Aerosol evolution 

Figure 4.3 depicts the predicted temperature profile as a function of the axial 
distance x and the jet radius r. The initial distribution is characterized by a hot, 
isothermal core and a somewhat colder bypass region. At the nozzle exit plane the 
initial values over the core and the bypass cross sections are assumed to be uniform. 
The inner bypass zone is heated up due to mixing of heat out of the core in the very 
near field. The plume temperature decays within the mixing region as the ambient 
air is entrained into the plume. 

Figure 4.4 presents the profile of soot emitted by the engine and transported as 
a passive tracer in the jet flow field. We assume that there is no soot present in 
the ambient air. The initial profile of the exhaust particulates decays, similar to the 
trend of the temperature, but a strong perturbation of the background atmosphere 
with respect to soot is observed, even at the simulated plume length of x = 1 km. 

Figure 4.5 shows the homogeneous nucleation rate along the plume axis. The 
present modeling results reveal that the typical germ radii are 0.3nm with 1.5 sulfuric 
acid and 3 water molecules per germ on average. The freshly nucleated acid aerosols 
have H 2 S0 4 weight fractions ranging from 0.7-0.8 depending on germ radius. The 
nucleation rate reaches its maximum by 100 m from the exit plane , and then it 
plummets as gaseous H2SO4 is depleted. 
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Figure 4.6 depicts the evolution of number densities of H 2 S0 4 -H 2 0 aerosols, soot 
particles, and soot coated with H 2 S0 4 -H 2 0 along the plume axis with and without 
coagulation effects. High levels of volatile H 2 S0 4 -H 2 0 droplets are predicted by 100 
m past the nozzle exit. After reaching its maximum, the amount of volatile droplets 
decreases due to both dilution and coagulation as manifested in the figure. The 
predicted number densities of bare soot and soot coated with H 2 S0 4 -H 2 0 coincide 
with those of mixed soot without coagulation effects. This suggests that we can 
combine the bare soot and those coated with H 2 S0 4 -H 2 0 into one parameter to 
save computing time. Here the generation path of H 2 S0 4 -H 2 0 coated soot particles 
through adsorption is neglected. 

Figure 4.7 shows the mixing ratio of gaseous sulfur species and H 2 S0 4 -H 2 0 aerosols. 
It clearly shows that the binary homogeneous nucleation effectively removes H 2 S0 4 
species from the gas phase. By the time of 1 second (c.f. Figure 3.10), the gas- phase 
H 2 S0 4 species is almost depleted. It indicates that any S0 2 that is oxidized to SO3 
and subsequently to sulfuric acid (H 2 S0 4 ) would form new small aerosols in the first 
second after emission. This estimate is comparable with the simulation results by 
Miake-Lye et al. [1994] and Karcher [1996], and is also consistent with the prediction 
of < 1.6 sec by Arnold et al.’s [1998] measurements. 
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4.5.2 Soot Activation 

Figure 4.8 shows the time evolution of the average surface coverage of sulfur 
species on soot particles along the plume axis. Prior to the onset of homogeneous bi- 
nary nucleation, adsorption dominates. This process can be viewed as prompt binary 
heterogeneous nucleation of H2SO4 and H2O on the surface. After homogeneous bi- 
nary nucleation, gaseous H2SO4 is depleted and the scavenging dominates the surface 
activation. 

4.5.3 Effects of sulfur emission from engine 

To explore the effect of SO* emissions at the exhaust exit on the formation, growth, 
and chemical processing of near field aerosols, parametric studies are performed. Since 
direct measurements of the level of S(VI) (SO3 and H 2 S0 4 ) at the exhaust exit have 
not been reported, the percentage conversion of fuel sulfur into SO3 in the engine 
was varied between 2% and 20%. This range encompasses both Brown et al’s [1996b] 
model estimates for the Concord and a smaller subsonic engine as well as experimental 
measurements of gas turbine exhausts [Hunter, 1982; Harris, 1990]. 

Figure 4.9 presents the maximum conversion efficiency of emitted SO* to H2SO4 
along the plume axis versus percentage conversion of fuel sulfur to SO3 in the engine. 
[SO3]o/[S0x]o is used as a free parameter since direct measurements of SO3 at the 
engine exit are not available. The subscript 0 indicates initial plume concentration 
at the nozzle exit plane. For low levels of SO3 emission with [S03]o/[S0x]o < 1%, 
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[i^SO^max/ISOxJo depends strongly on [OH] 0 , because SO 3 levels produced via S0 2 
+ OH are greater than or at least similar to [SOaJo- Limited by the S0 2 4- OH 
reaction rate, [# 2 ‘S 04 ] m ax/[‘SO I ]o is bounded between 0 . 1 % and 2 % for [OH ] 0 levels 
between O.lppmv and lOppmv. Direct emission of H 2 SO 4 could enhance this limit. 
However, when [50 3 ] 0 /[50 x ]o > 1%, [/^S^max/fSOxJo increases dramatically be- 
cause the fast SO 3 + H 2 O reaction produces H 2 S0 4 directly, bypassing the slower 
S0 2 + OH reaction pathway and eliminating the [OH]o dependence. Direct emission 
of H 2 S0 4 and possible heterogeneous losses of SO 3 and H 2 SO 4 due to adsorption on 
soot particles were not considered here. The computed conversion range is consistent 
with the estimated value of > 0.4% based on measurements by Frenzel and Arnold 
[1994], 

Figure 4.10 shows the accumulative number densities (as would be detected by 
a condensation nucleus counter) of volatile H 2 S0 4 -H 2 0 aerosols with radius greater 
than 3.5 nm along the jet plume axis. Profiles are shown assuming that 2%, 5%, 10%, 
and 20 % of the fuel sulfur is converted to SO 3 in the engine. The predicted cumulative 
number density of volatile aerosols is shown to increase with the S0 3 emission. As 
more SO 3 is emitted from the nozzle exit, more gas-phase H 2 S0 4 is produced as shown 
in Figure 4.9. This leads to higher peak homogeneous binary nucleation rates and 
therefore more volatile H 2 SO 4 -H 2 O embryos are generated in the early jet regime. 
Elevated levels of gaseous H 2 SO 4 also promote the condensation process which causes 
volatile acid droplets to grow into larger size. This is reflected in the earlier and 
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sharper rise of the average aerosol radius for higher SO3 emission as depicted in Figure 
4.11. The average radius of aerosols is defined as r VAVe = • 47r), where 

A Vt tot and N Vjto t are the total surface area and number density of volatile particles in 
a grid, respectively. When all the SOx emitted from nozzle are in the form of SO2, 
the cumulative number densities are too small to be shown in Figure 4.10. In the late 
stage of the jet plume the evolution of number densities are shown to be dominated 
by dilution. 

In Figure 4.12, the evolution of relative humidity (RH) along the plume axis is 
illustrated. The relative humidity first increases due to rapid temperature decrease 
caused by the entrainment of cold ambient air. The water vapor reaches supersatura- 
tion at about 100 m after exit. The level of water vapor supersaturation achieves its 
maximum by 250 m, then it falls off due to condensation of water vapor onto aerosols 
until RH as 100%. Hereafter, particles begin to evaporate water due to entrainment. 
The net effect is that RH decreases slightly. With higher conversion of fuel sulfur to 
oxidized sulfur in the engine, more and larger volatile aerosols are generated leading 
to higher levels of surface area density to be presented in Figure 4.13. Hence en- 
hanced condensation due to larger surface area available causes RH to decrease to 
100% sooner as depicted in Figure 4.12. 

A natural and important consequence of the generation of high number density of 
volatile H2SO4-H2O aerosols is a significant increase in aerosol surface area density. 
Figure 4.13 presents the volatile aerosol surface area density along the plume axis 



78 


for the conditions with 0, 2 %, 5%, 10%, and 20% emission of SO* as S0 3 at the 
exhaust exit. Homogeneous binary nucleation first causes a marked increase of surface 
area density of volatile acid aerosols. After reaching its maximum at about 50 m, 
nucleation is slowed down, and surface area density is decreased by coagulation. Once 
the plume becomes supersaturated with respect to water, the surface area increases 
due to rapid water condensation on droplets with radii larger than nanometer size. 
After the plume relative humidity decreases to about 100 %, condensation ceases and 
the surface area density decreases which is dominated by plume dilution and aerosol 
coagulation. When S0 3 emission is zero at the nozzle exit, the limited amount of 
H2SO4 is depleted by adsorption and homogeneous binary nucleation before the plume 
becomes saturated. The sizes of volatile aerosols stay small and water uptake will not 
be important due to the Kelvin barrier. Hence the surface area does not exhibit a 
second increase period which is caused by water condensation. For higher percentage 
emission of SOx as SO3, more H2SO4 is generated, which leads to higher homogeneous 
binary nucleation rate and thus generates more volatile acid aerosols. This results 
in larger increase in surface area early in the plume. Higher levels of gas phase 
H2SO4 also enhance H2SO4 condensation and accordingly induce more H2O uptake, 
which makes particles grow into larger size as manifested in Figure 4.11. Therefore 
a second increase period becomes more profound with higher levels of SO3 emission 
from the nozzle exit. The computed surface area density levels are in the range of 
0(10 3 )-0(10 5 ) /im 2 cm' 3 which are consistent with the predictions by Karcher [1996]. 
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As shown in Figure 4.14, the level of chemical activation of soot particles increases 
as fuel sulfur conversion to S 0 3 in the engine increases. When the exhaust from the 
nozzle exit contains no SO3, 6 is about 0.3 at the end of jet plume regime. If 2% of 
fuel sulfur is converted to SO3 in the engine, 6 would get to nearly 0.8 at the end of 
jet. When the conversion fraction is > 5%, 6 reaches nearly 1. 

In the near-field jet plume, the computed results indicate that only l%-5% of 
exhaust SO3 is adsorbed on soot during the time required for chemical buildup of 
maximum H2SO4 (10 ms) in the gas phase for the range of SO3 emission studied 
here. Hence the adsorption of SO3 is not a limit to its gas-phase chemical conversion 
to H2SO4. The H2SO4 generated in the gas phase is adsorbed by soot until the 
depletion of available H2SO4 by homogeneous binary nucleation and condensation 
around a plume age of 0.1 sec. The model results show that approximately 3%-25% 
of the chemically produced H2SO4 is adsorbed by soots. Hence we can conclude that 
H2SO4 is the major species adsorbed on soot particles compared to the adsorption 
of SO3 on soot. However, the amount of H2SO4 adsorbed on soot particles is still 
relatively small and it will not limit the H2SO4-H2O binary homogeneous nucleation 
and H2SO4 condensation. 

Based on simulation of combustion chemistry, Brown et al. [1996b] concluded that 
the conversion fraction of fuel sulfur to S(VI) (SO3+H2SO4) does not scale linearly 
with fuel sulfur content (FSC). The conversion is lower for higher FSC. 
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4.5.4 Threshold temperature for contrail formation 

Jet-Al fuel is a standard aviation kerosene with hydrogen mass fraction of m//= 0.14 
and specific combustion heat of Q = 43 x 10 6 J/kg [Schumann, 1996]. For a B747 air- 
craft burning kerosene with a thrust F = 31.1 x 10 3 N per engine [Schumann, 1996], 
a fuel flow rate m = 0.785 kg/s, and the air velocity V = 237 m/s, the propulsion effi- 
ciency and the H 2 O emission index are estimated to be 77 = 0.216 and EIh 2 o — 1-25, 
respectively. With the ambient pressure P a = 0.2361 atm and relative humidity 
RH a = 30%, we compute the threshold ambient temperature Tic to be 221.51 K. 
For a higher ambient relative humidity RH a = 50%, Tic= 222.48 K. As expected, 
ambient air with higher relative humidity facilitates contrail formation. For r? = 0 
and RH = 50%, Tic is 220.19 K which is about 2 K smaller than the situation with 
•q = 0.216. Hence, the Appleman’s criterion would give lower threshold temperatures 
since it did not consider the conversion of combustion heat into kinetic energy in the 
wake vortex regime. 

Given the overall propulsion efficiency of an aircraft and the type of fuel burnt, the 
threshold temperature turns out to be a function of ambient relative humidity and 
pressure (altitude), as seen from equations (4.65), (4.66), and (4.68). In Figure 4.15, 
the computed threshold temperature for contrail formation is plotted versus ambient 
pressure for ambient relative humidities of 0, 50, and 100%. The propulsion efficiency 
r} is specified to be 0.3, which is the approximate value for modern commercial jet 
engines. It clearly indicates that contrails tend to form at higher ambient pressure or 
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higher relative humidity. The threshold temperature is seen sensitive to the aircraft 
propulsion efficiency, and contrails tend to form at a higher propulsion efficiency. 
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temperature (K) 

Figure 4.2: Schematic showing threshold condition for contrail formation. It shows the 
curve of the liquid water saturation pressure P sa t,Uq{T) and ice saturation pressure 
Ptat,ice(T ) versus temperature and the mixing lines. T L c is the threshold ambient 
temperature for contrail formation. 
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Figure 4.3: Spatial distribution of temperature in the flow field of a B747 jet engine. 
The radial distance, r, and the downstream distance, x, are in unit m. 
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Figure 4.4: Spatial distribution of the total number density of soot particles in the 
flow field of a B747 jet engine. The radial distance, r, and the downstream distance, 
x, are in unit m. 



number density (#/cm 3 ) 


87 
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Figure 4.6: Number densities of volatile H 2 S0 4 -H 2 0 droplets, soot particles, and soot 
coated with H 2 SC>4-H 2 0 versus downstream distance with and without coagulation 
effects. 
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Figure 4.7: Mixing ratio of gaseous sulfur species and sulfate aerosols versus down- 
stream distance. The homogeneous nucleation rate Jhom in unit of ppmv/sec is also 
shown. Dashed line represents the evolution of gaseous H 2 SO 4 species when nucle- 
ation is included as compared to the pure dilution (solid line). 
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Figure 4.8: Average surface coverage of oxidized sulfur species for a soot versus plume 
age. It is composed of coverage by adsorption of gaseous H 2 SO 4 and by scavenging 
of volatile H2SO4-H2O aerosols. 
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Figure 4.10: Cumulative number densities of volatile H2SO4-H2O aerosols with radius 
greater than 3.5 nm along the jet plume axis, assuming that 2%, 5%, 10%, and 20% 
of the fuel sulfur is converted to SO3 in the engine. 
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Figure 4.12: The evolution of relative humidity along the plume axis. Profiles are 
shown assuming 0, 2 %, 5 %, 10 %, and 20% of the fuel sulfur is converted to SO3 in 
the engine. 
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Figure 4.13: Volatile H 2 S0 4 -H 2 0 aerosol total surface area density along the jet 
plume axis. Profiles are shown assuming 0, 2%, 5%, 10%, and 20% of the fuel sulfur 
is converted to SO 3 in the engine. 
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ambient pressure (atm) 

Figure 4.15: Threshold temperature for contrail formation versus ambient pressure 
for ambient relative humidities (RH) of 0, 50, and 100%. Contrails are predicted to 
form when the ambient temperature is below the curve. Threshold temperatures are 
also shown for an engine propulsion efficiency of 0.2 and 0.4. 
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Chapter 5 


Conclusions 


A detailed chemical mechanism between aircraft exhaust and ambient air in near- 
field jet plumes is compiled. The dominant oxidation kinetic pathways are identified 
through sensitivity analysis of major species with respect to gas phase reaction steps. 
The spatial distributions of major chemical species in the near-field jet plume are com- 
puted and analyzed. The potential impact of subgrid turbulent mixing on chemistry 
is investigated and found to be of negligible importance. 

In the vortex regime of the wake, large eddy simulations (LES) may be required 
to track the highly three-dimensional fluid dynamics. The expensive computation 
on fluid dynamics makes it impossible to incorporate detailed chemistry. A reduced 
mechanism is therefore required. An eight-step reduced mechanism is developed which 
captures the major chemical pathways in the near-field jet plume. Results obtained 
with the reduced chemistry show good agreement with those predicted by the detailed 
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chemistry. However, using the eight-step reduced chemistry has cut down computer 
CPU time by a factor of more than 3.5 for the near-field simulation. 

A comprehensive aerosol model is developed to simulate microphysical processes 
including sulfuric acid and water binary homogeneous nucleation, aerosol coagulation, 
non-equilibrium heteromolecular condensation, and sulfur-induced soot activation as 
water condensation nuclei. These microphysical processes are coupled with gas- phase 
chemistry and axisymmetric turbulent jet flow. 

The modeling results indicate that a large number of subnanometer-sized volatile 
H2SO4-H2O embryos are generated in the near-field plume. These embryos grow 
in size by condensation in the early jet regime and by self coagulation later. Soot 
particles can be activated into water condensation nuclei in the jet plume for possible 
contrail formation through both heterogeneous nucleation and scavenging of volatile 
H2SO4-H2O aerosols. 

The gas phase conversion of SO2 to H2SO4 is controlled by OH driven oxidation 
and it is limited to less than 1%. However, the amount of H2SO4 produced will in- 
crease via fast reaction of SO3 with H2O directly if SO3 emission due to fuel sulfur 
oxidation to SO3 in the engine is considered. Higher levels of H2SO4 in turn result in 
more sulfate aerosol embryos formed in the early jet stage, and also enhance conden- 
sation of H2SO4 and H2O onto the sulfate aerosols in this regime. Consequently, S 0 3 
emission from the engine nozzle exit increases the cumulative number density and 
surface area density of volatile sulfate aerosols. SO3 emission also promotes chemical 
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activation of soot particles with enhanced acid coating. This suggests the need for 
measurements of the partitioning of sulfur species (S0 2 and S0 3 ) at the exhaust exit. 

The threshold condition for ice contrail formation behind an aircraft depends on 
the ambient pressure, the ambient relative humidity, the overall propulsion efficiency 
of the aircraft, and the fuel properties such as the fuel hydrogen mass fraction and 
the specific combustion heat. Higher values of ambient relative humidity or pressure 
facilitate contrail formation. The threshold temperature is found sensitive to the 
overall propulsion efficiency of the aircraft. 

The prediction of particle and species distribution in the near-field plume can 
serve as input for the follow-on larger-scale model. The present work also provides an 
estimate of the aerosol surface area density for subsequent evaluation of heterogeneous 


oxidation mechanisms in the wake. 
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